Building an Exome Analysis Pipeline From FASTQ to Annotated Variants

By the end of this guide you will have taken paired-end exome FASTQ files and produced three artifacts: a coordinate-sorted CRAM aligned to GRCh38 with duplicates marked and base qualities recalibrated, a filtered VCF of single-nucleotide variants and short insertions and deletions restricted to the capture target plus flanking bases, and an annotated table with gene, consequence, population frequency, and in-silico prediction scores per variant. You will also have a coverage report and a concordance report against a reference truth set, which is the part most tutorials skip and the part that tells you whether the rest is trustworthy. You need a Linux machine with 8 or more cores and at least 32 GB of RAM, and roughly 200 GB of free disk per sample if you keep intermediates. You need the following tools on it.

• bwa-mem2 (or bwa)
• samtools 1.18+
• GATK 4.5+
• bcftools
• fastp
• mosdepth
• Ensembl VEP 110+
• hap.py or rtg vcfeval for benchmarking

You also need the BED file describing the capture kit your library was built with. If you do not know which kit was used, get it from whoever ran the library prep before you do anything else, because target definition affects every downstream number.

Exome sequencing captures the roughly 1 to 2 percent of the genome that is protein-coding, which is where a large share of known Mendelian disease alleles sit, and does so at a fraction of the cost and data volume of whole genomes ¹. The approach was demonstrated as a diagnostic tool early, with exome capture plus massively parallel sequencing identifying a causal variant in a single patient with a suspected Mendelian condition ². That efficiency comes with a specific cost: capture is uneven, off-target reads are wasted, and the boundaries of your target define what you can and cannot say.

1. Inspect the reads and trim only what needs trimming

Start by looking at the raw data rather than assuming it is fine. Run fastp, which gives you per-base quality, adapter content, duplication estimate, and insert size distribution in one pass, and writes a JSON you can parse later.

fastp \
  -i sample_R1.fastq.gz -I sample_R2.fastq.gz \
  -o sample.trim_R1.fastq.gz -O sample.trim_R2.fastq.gz \
  --detect_adapter_for_pe \
  --cut_tail --cut_tail_mean_quality 20 \
  --length_required 50 \
  --thread 8 \
  --json sample.fastp.json --html sample.fastp.html

We deliberately avoid aggressive quality trimming. BWA-MEM soft-clips poor-quality ends on its own, and GATK’s base quality score recalibration models systematic error better than a sliding-window trimmer does. What matters is adapter removal. Adapter sequence left on short-insert fragments produces spurious soft clips and, in the worst case, false indel calls at read ends. Check the insert size histogram in the HTML report: a median insert below about 150 bp on 2x150 chemistry means substantial read-through into adapter and overlapping mates, and you should expect inflated duplicate rates.

2. Align to GRCh38 with alt-aware BWA-MEM

Use GRCh38 with the analysis set that includes decoy and alt contigs. The decoys absorb reads from sequences absent in the primary assembly, which otherwise pile up in paralogous regions and generate false heterozygous calls.

bwa-mem2 mem -t 16 -K 100000000 -Y \
  -R '@RG\tID:S1.L1\tSM:S1\tLB:lib1\tPL:ILLUMINA\tPU:FLOWCELL.1.S1' \
  GRCh38_full_analysis_set_plus_decoy_hla.fa \
  sample.trim_R1.fastq.gz sample.trim_R2.fastq.gz \
| samtools fixmate -m -u - - \
| samtools sort -@ 4 -m 2G -o sample.sorted.bam -
samtools index -@ 8 sample.sorted.bam

The read group line is not decoration. GATK refuses to run without it. The PU (platform unit) field is what BQSR uses to model error per lane. For that reason, if you sequenced across multiple lanes, align each lane separately with its own read group and merge. The -Y flag uses soft clipping for supplementary alignments, which keeps sequence available for structural inspection later. The -K 100000000 setting fixes the batch size so output is deterministic regardless of thread count, which matters when you want to reproduce a result exactly.

If you are using alt contigs, run the bwa-postalt.js script from the bwakit distribution before sorting. Skipping it leaves reads with low mapping quality in regions like HLA and some immunoglobulin loci, and you will see coverage holes there that look like capture failures but are mapping artifacts.

3. Mark duplicates and recalibrate base qualities

Exome libraries have higher duplicate rates than genomes because capture reduces library complexity. Mark rather than remove, so you can revisit the decision.

gatk MarkDuplicatesSpark \
  -I sample.sorted.bam \
  -O sample.dedup.bam \
  -M sample.dup_metrics.txt \
  --optical-duplicate-pixel-distance 2500 \
  --conf 'spark.executor.cores=16'

gatk BaseRecalibrator \
  -I sample.dedup.bam -R GRCh38_full_analysis_set_plus_decoy_hla.fa \
  --known-sites Homo_sapiens_assembly38.dbsnp138.vcf \
  --known-sites Mills_and_1000G_gold_standard.indels.hg38.vcf.gz \
  --known-sites Homo_sapiens_assembly38.known_indels.vcf.gz \
  -L targets.padded100.bed \
  -O sample.recal.table

gatk ApplyBQSR \
  -I sample.dedup.bam -R GRCh38_full_analysis_set_plus_decoy_hla.fa \
  --bqsr-recal-file sample.recal.table \
  -O sample.recal.cram

Two notes on the flags. The optical duplicate pixel distance of 2500 is correct for patterned flow cells (NovaSeq, HiSeq X). The default of 100 is for unpatterned instruments and will undercount duplicates on a NovaSeq. And restricting BQSR training with -L to the padded target keeps the model from being dominated by off-target reads, which have a different error profile.

Build the padded BED once and reuse it everywhere:

bedtools slop -i targets.bed -g GRCh38.genome -b 100 \
| bedtools sort | bedtools merge > targets.padded100.bed

We use 100 bp of padding rather than the more common 50 or 10. Capture probes pull in fragments that extend beyond the probe footprint, so coverage falls off gradually, and canonical splice acceptor and donor sites sit just outside many coding exon definitions. Padding wider costs you a slightly larger callset to filter and buys you the splice region.

4. Check coverage before you call anything

Coverage is where exomes fail, and knowing the failure map up front saves you from over-interpreting a negative result. Different capture platforms differ substantially in what they target and how evenly they cover it, with meaningful disagreement in the variants they detect in the same sample ³.

mosdepth --by targets.padded100.bed --fast-mode \
  --thresholds 1,10,20,30 -t 4 sample sample.recal.cram

Read two numbers from sample.thresholds.bed.gz and sample.mosdepth.summary.txt: the fraction of target bases at 20x or more, and the mean target depth. For germline single-nucleotide variant calling in a diploid sample, we want 90 percent or more of target bases at 20x and a mean of 80x to 100x. Mean depth alone is misleading because GC-rich first exons are systematically under-covered, so a sample with a 120x mean can still have thousands of bases below 10x. Extract the low-coverage intervals and keep them alongside your VCF:

zcat sample.thresholds.bed.gz \
| awk 'NR>1 && $5/($3-$2) < 0.9' \
| cut -f1-4 > sample.lowcov_regions.bed

Any gene on that list is a gene where the absence of a variant call means nothing.

5. Call variants with HaplotypeCaller in GVCF mode

Call in GVCF mode even for a single sample. A GVCF records confidence at reference positions as well as variant sites, which lets you distinguish “reference” from “no data” and lets you joint-genotype later if you add family members.

gatk --java-options "-Xmx16g" HaplotypeCaller \
  -R GRCh38_full_analysis_set_plus_decoy_hla.fa \
  -I sample.recal.cram \
  -L targets.padded100.bed \
  --interval-padding 0 \
  -ERC GVCF \
  -G StandardAnnotation -G AS_StandardAnnotation \
  -O sample.g.vcf.gz

gatk GenotypeGVCFs \
  -R GRCh38_full_analysis_set_plus_decoy_hla.fa \
  -V sample.g.vcf.gz \
  -L targets.padded100.bed \
  -O sample.raw.vcf.gz

HaplotypeCaller performs local de novo assembly of haplotypes in active regions, which is why it outperforms position-based callers on indels and in regions with clustered variation. Request allele-specific annotations (AS_StandardAnnotation) now even if you filter on the site-level ones, because you cannot recover them later without recalling.

One structural point worth internalizing: the choice of caller and aligner changes your answer. Comparisons of independent exome pipelines run on the same reads have found concordance well below what most users assume, with a substantial fraction of single-nucleotide variant calls unique to one pipeline and indel agreement considerably worse ⁴. Later assessments using both real exome data and simulated reads with known truth reached the same conclusion, that sensitivity and precision vary by caller and are worst for indels ⁵. Treat any single pipeline’s output as one estimate. It is not ground truth.

6. Filter

With one sample and an exome, you do not have enough variants for GATK’s Variant Quality Score Recalibration to train reliably on indels. Use hard filters, applied separately to SNVs and indels because they have different quality distributions.

gatk SelectVariants -V sample.raw.vcf.gz -select-type SNP -O snps.vcf.gz
gatk SelectVariants -V sample.raw.vcf.gz -select-type INDEL -select-type MIXED -O indels.vcf.gz

gatk VariantFiltration -V snps.vcf.gz -O snps.filt.vcf.gz \
  --filter-name QD2 --filter-expression "QD < 2.0" \
  --filter-name QUAL30 --filter-expression "QUAL < 30.0" \
  --filter-name SOR3 --filter-expression "SOR > 3.0" \
  --filter-name FS60 --filter-expression "FS > 60.0" \
  --filter-name MQ40 --filter-expression "MQ < 40.0" \
  --filter-name MQRS-12.5 --filter-expression "MQRankSum < -12.5" \
  --filter-name RPRS-8 --filter-expression "ReadPosRankSum < -8.0"

gatk VariantFiltration -V indels.vcf.gz -O indels.filt.vcf.gz \
  --filter-name QD2 --filter-expression "QD < 2.0" \
  --filter-name QUAL30 --filter-expression "QUAL < 30.0" \
  --filter-name FS200 --filter-expression "FS > 200.0" \
  --filter-name RPRS-20 --filter-expression "ReadPosRankSum < -20.0"

gatk MergeVcfs -I snps.filt.vcf.gz -I indels.filt.vcf.gz -O sample.filt.vcf.gz

Then add a genotype-level filter for depth and allele balance, which catches a category of error the site-level filters miss:

bcftools filter -S . -e 'FMT/DP<10 | FMT/GQ<20' sample.filt.vcf.gz \
| bcftools view -f PASS -e 'GT="mis"' -Oz -o sample.pass.vcf.gz
bcftools index sample.pass.vcf.gz

If you have three or more exomes sequenced on the same kit and instrument, switch to VQSR for SNVs, which uses the known-sites resources to learn the boundary between true and false calls rather than imposing fixed thresholds. Keep hard filters for indels regardless.

7. Annotate with VEP

Annotation determines what you can interpret, so pin the versions of both the tool and the cache and record them.

vep --cache --offline --dir_cache /refs/vep --assembly GRCh38 \
  --species homo_sapiens --cache_version 110 \
  --input_file sample.pass.vcf.gz --format vcf \
  --output_file sample.vep.vcf.gz --vcf --compress_output bgzip \
  --everything \
  --pick_order rank,canonical,tsl,biotype,length --pick \
  --plugin CADD,/refs/cadd/whole_genome_SNVs.tsv.gz \
  --plugin SpliceAI,snv=/refs/spliceai/snv.vcf.gz,indel=/refs/spliceai/indel.vcf.gz \
  --custom /refs/gnomad/gnomad.exomes.v4.1.sites.vcf.bgz,gnomADe,vcf,exact,0,AF,AF_popmax,nhomalt \
  --fork 8 --buffer_size 20000 --stats_file sample.vep.html

Use --pick for a per-variant summary but keep the full transcript-level output too, because the highest-ranked consequence on an arbitrary transcript can hide a more relevant consequence on the transcript that matters for a given gene. The --custom gnomAD annotation is what makes filtering tractable: population allele frequency is the single most effective filter for narrowing thousands of coding variants to a reviewable set, and standard disease-gene discovery strategies rest on it together with consequence class and inheritance model ⁶. A practical first pass is to keep variants with gnomAD allele frequency below 0.001 and a consequence class of missense or worse, then layer in zygosity and gene lists.

8. Benchmark against a truth set

Run a reference sample through the exact pipeline you just built. HG002 or NA12878 from Genome in a Bottle, with the corresponding high-confidence VCF and BED, intersected with your capture target.

bedtools intersect -a HG002_GIAB_highconf.bed -b targets.padded100.bed > eval.bed

hap.py HG002_GIAB_highconf.vcf.gz hg002.pass.vcf.gz \
  -f eval.bed -r GRCh38_full_analysis_set_plus_decoy_hla.fa \
  --engine vcfeval --threads 8 -o hap_hg002

Read precision and recall separately for SNVs and indels from hap_hg002.summary.csv. On a well-behaved exome at 100x we expect SNV recall and precision above 0.99 inside high-confidence regions, and indel recall in the 0.92 to 0.97 range with lower precision. If your indel numbers are far below that, suspect adapter contamination, over-trimming, or a missing fixmate step. This benchmarking loop is the only way to know what your specific pipeline does, and the value of running it is exactly why automated single-command exome pipelines were built with multiple interchangeable callers and a concordance step ⁷.

Common problems

Coverage holes in GC-rich first exons are the most common cause of a false negative, and they are a property of hybridization capture rather than a bug in your code. First exons are often more than 65 percent GC, probes there hybridize poorly, and PCR amplification during library prep compounds the bias. Always read a negative result against your low-coverage BED before concluding a gene is clear.

Off-target reads consume sequencing capacity you paid for. A typical figure of merit is the fraction of bases falling inside the padded target, and 60 to 80 percent is normal depending on kit. Below about 50 percent, suspect a mismatch between the BED file and the kit used. Other causes are insufficient hybridization time or a library with high adapter dimer content.

Indel calls are systematically less reliable than SNV calls, especially in homopolymers and short tandem repeats. Left-align and normalize before comparing anything (bcftools norm -f ref.fa -m -both), because two callers can represent the same indel at different coordinates and appear discordant when they agree. Where an indel matters, look at the reads in IGV rather than trusting the quality score.

Copy number variants are largely invisible to this pipeline. Exome capture gives you depth in disjoint intervals, so single-exon deletions and duplications require a dedicated read-depth CNV caller with a panel of normals sequenced on the same kit, and the results are noisy in small events. Similarly, variants in deep introns, promoters, and untranslated regions are outside the target, and structural rearrangements with breakpoints in intronic sequence will be missed. These blind spots are inherent to the method and have been described since exome sequencing moved into disease studies ⁸.

Pipeline drift is a quiet failure mode. Reference builds, VEP caches, gnomAD versions, and GATK releases all change behavior, and a variant that passed last year can fail today. Record tool versions and reference checksums in the output, containerize the whole thing, and rerun your benchmark sample after any upgrade. Published comparisons of exome pipelines consistently attribute a large share of disagreement to software and parameter choices rather than to the data, which is why version pinning and re-benchmarking matter more than picking the nominally best tool ⁹.

Finally, interpretation is a separate discipline from computation. Producing a short list of rare, predicted-damaging variants in genes plausibly related to a phenotype is a bioinformatics result. The list is not a clinical finding. Confirming a candidate, judging its pathogenicity against curated evidence, and deciding what it means for a person requires a clinical geneticist or genetic counselor working with a clinically validated laboratory report. A research-grade pipeline like the one above should not be used to make health decisions on its own. A published protocol walking through the same steps end to end is a useful cross-check on your implementation ¹⁰.

Footnotes

1. J. K. Teer, J. C. Mullikin. Exome sequencing: the sweet spot before whole genomes. Human Molecular Genetics, 2010. https://doi.org/10.1093/hmg/ddq333 ↩

2. Murim Choi, Ute I. Scholl, Weizhen Ji, et al. Genetic diagnosis by whole exome capture and massively parallel DNA sequencing. Proceedings of the National Academy of Sciences, 2009. https://doi.org/10.1073/pnas.0910672106 ↩

3. 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 ↩

4. Jason O’Rawe, Tao Jiang, Guangqing Sun, et al. Low concordance of multiple variant-calling pipelines: practical implications for exome and genome sequencing. Genome Medicine, 2013. https://doi.org/10.1186/gm432 ↩

5. Manojkumar Kumaran, Umadevi Subramanian, Bharanidharan Devarajan. Performance assessment of variant calling pipelines using human whole exome sequencing and simulated data. BMC Bioinformatics, 2019. https://doi.org/10.1186/s12859-019-2928-9 ↩

6. PN Robinson, P Krawitz, S Mundlos. Strategies for exome and genome sequence data analysis in disease‐gene discovery projects. Clinical Genetics, 2011. https://doi.org/10.1111/j.1399-0004.2011.01713.x ↩

7. Yunfei Guo, Xiaolei Ding, Yufeng Shen, et al. SeqMule: automated pipeline for analysis of human exome/genome sequencing data. Scientific Reports, 2015. https://doi.org/10.1038/srep14283 ↩

8. Zuoheng Wang, Xiangtao Liu, Bao-Zhu Yang, et al. The Role and Challenges of Exome Sequencing in Studies of Human Diseases. Frontiers in Genetics, 2013. https://doi.org/10.3389/fgene.2013.00160 ↩

9. Zeeshan Ahmed, Eduard Gibert Renart, Saman Zeeshan. Genomics pipelines to investigate susceptibility in whole genome and exome sequenced data for variant discovery, annotation, prediction and genotyping. PeerJ, 2021. https://doi.org/10.7717/peerj.11724 ↩

10. Narendra Meena, Praveen Mathur, Krishna Medicherla, et al. A Bioinformatics Pipeline for Whole Exome Sequencing: Overview of the Processing and Steps from Raw Data to Downstream Analysis. BIO-PROTOCOL, 2018. https://doi.org/10.21769/bioprotoc.2805 ↩