How to Call Variants From Your Own Whole-Genome Sequencing Data

By the end of this you will have a per-sample small-variant VCF (roughly 4.5 to 5 million SNVs and indels against GRCh38), a structural variant call set, a repeat-expansion genotype table, an annotated TSV you can query, and a set of quality metrics that tell you whether any of it is trustworthy. You need: your raw data (FASTQ.gz or an aligned CRAM/BAM, 30x PCR-free Illumina is the assumption throughout), a machine with 32 cores and 64 GB of RAM, about 1 TB of scratch disk, Docker or Singularity, and a day of wall-clock time if you are running end to end from FASTQ. Everything below is command-line. Nothing here is a diagnosis, and nothing here substitutes for a clinical laboratory report.

1. Inventory and verify what you received

Before any compute, confirm the files are intact and know which flavor you have.

md5sum -c checksums.md5
samtools quickcheck -v *.cram && echo "cram ok"
samtools view -H sample.cram | grep '^@RG'
samtools view -H sample.cram | grep -E '^@SQ' | head -3

The @RG line tells you the read group, platform, and library. If PL is not ILLUMINA or PL:ONT/PACBIO, find out before you pick a caller. The @SQ lines tell you the reference the file was aligned to. Watch for the SN:chr1 versus SN:1 distinction and for the presence of chrUn_* decoys and *_alt contigs, because that determines which reference FASTA you must use downstream. A CRAM without the exact matching reference is unreadable. It is not approximately readable.

Coverage first, because low coverage invalidates everything after it:

mosdepth --by 1000 --fast-mode --threads 4 sample sample.cram
zcat sample.regions.bed.gz | awk '{s+=$4*($3-$2); n+=($3-$2)} END {print s/n}'

For germline small-variant calling on short reads, we want mean autosomal coverage at or above 30x and at least 90 percent of the callable genome at 20x or better. Sensitivity for heterozygous sites falls off sharply below 20x because you need both alleles sampled several times. Clinical sequencing guidance treats depth per site, not mean depth, as the number that matters ¹.

2. Fix the reference genome before you fix anything else

Use the GRCh38 no-ALT analysis set with decoy sequences (GCA_000001405.15_GRCh38_no_alt_analysis_set.fna). The reason is specific: ALT contigs create multi-mapping reads that depress mapping quality across the MHC and other segmental duplications, and most callers silently drop MAPQ 0 reads. The decoy sequences absorb reads from sequences absent in the primary assembly (EBV, unplaced repeats) that would otherwise pile onto false homologous loci.

If your data arrived aligned to hs37d5 or GRCh37, realign from the CRAM rather than lifting the VCF over. Liftover of coordinates is lossy at exactly the loci you care about, and reference-allele flips silently invert genotypes. Realignment costs a few hours of CPU. Liftover costs you correctness you cannot audit.

T2T-CHM13 is a better assembly and resolves regions GRCh38 cannot. Most annotation resources (gnomAD frequencies, ClinVar coordinates, VEP caches) are nonetheless keyed to GRCh38. We run GRCh38 as primary and CHM13 as a secondary call set when something interesting lands in a segmental duplication.

3. Align, if you are starting from FASTQ

bwa-mem2 mem -t 32 -K 100000000 -Y \
  -R '@RG\tID:L001\tSM:SAMPLE\tLB:lib1\tPL:ILLUMINA\tPU:FLOWCELL.1' \
  GRCh38_no_alt.fna R1.fastq.gz R2.fastq.gz \
| samtools sort -@ 8 -m 2G -T /scratch/tmp -o sample.sorted.bam -

gatk MarkDuplicates \
  -I sample.sorted.bam -O sample.md.bam -M sample.dup_metrics.txt \
  --OPTICAL_DUPLICATE_PIXEL_DISTANCE 2500

samtools view -@ 8 -C -T GRCh38_no_alt.fna -o sample.cram sample.md.bam
samtools index sample.cram

Two flags earn their keep. -K 100000000 fixes the batch size so output is deterministic regardless of thread count, which matters when you want to reproduce a call set exactly. -Y uses soft clipping for supplementary alignments, which structural variant callers need to find split reads.

Skip base quality score recalibration if you are running DeepVariant. BQSR was designed for GATK’s Bayesian model, and the neural-network callers learn base quality behavior from the data. If you are running GATK HaplotypeCaller, run BQSR. The functional-equivalence work behind the large sequencing consortia specified exactly which pipeline steps can differ between implementations while producing concordant calls, and the alignment and duplicate-marking steps are where divergence creeps in first ².

4. Call small variants

Our recommendation for 30x PCR-free Illumina germline data is DeepVariant. It converts pileups to images and classifies genotypes with a CNN, which means it does not need per-cohort hard-filter tuning, and it handles indels in homopolymer runs better than a hand-tuned hard filter.

docker run -v "$PWD":/data google/deepvariant:1.6.1 \
  /opt/deepvariant/bin/run_deepvariant \
  --model_type=WGS \
  --ref=/data/GRCh38_no_alt.fna \
  --reads=/data/sample.cram \
  --output_vcf=/data/sample.dv.vcf.gz \
  --output_gvcf=/data/sample.dv.g.vcf.gz \
  --num_shards=32 \
  --intermediate_results_dir=/data/tmp_dv

Roughly 2 to 4 hours on 32 cores at 30x. Always write the gVCF. You will want it later if you ever join your sample with family members through GLnexus, and the gVCF is the only record of which reference-matching positions were confidently called versus merely uncovered.

The alternatives, and when we would use them:

• GATK HaplotypeCaller in -ERC GVCF mode, then GenotypeGVCFs plus VQSR or hard filters. Use it if you need to match an existing consortium pipeline, since the functional-equivalence specification was written around this family of pipelines ².
• Sentieon’s DNAscope/DNAseq. Same mathematical models as GATK, reimplemented for speed, with identical or near-identical output and order-of-magnitude runtime reductions ³. Worth it if you are calling many samples on your own hardware.
• Strelka2 if you want something fast and CPU-light with good indel behavior and no GPU.

Set --model_type=WES for exome data, and use the PacBio or ONT models for long reads. Running the WGS model on exome data inflates false positives at capture boundaries where coverage cliffs are steep.

5. Normalize and filter

Raw caller output is not queryable. Normalize first, always, before comparing to anything:

bcftools norm -m -any -f GRCh38_no_alt.fna -Ou sample.dv.vcf.gz \
| bcftools norm -d exact -Oz -o sample.norm.vcf.gz
bcftools index -t sample.norm.vcf.gz

-m -any splits multiallelic records into biallelic rows. -f left-aligns indels against the reference. Without left-alignment, the same indel written two ways will not match ClinVar, gnomAD, or your friend’s VCF.

Then apply a conservative quality filter:

bcftools view -f PASS -i 'FMT/GQ>=20 && FMT/DP>=10' \
  -Oz -o sample.pass.vcf.gz sample.norm.vcf.gz

GQ 20 corresponds to a 1 percent posterior probability of the wrong genotype. Raising it to 30 buys precision and costs recall at heterozygous sites near the coverage floor. Which way you want that tradeoff depends on whether you are screening broadly or following up a specific locus.

Exclude the hard regions from any conclusion you plan to act on. Download the GIAB genome stratification BEDs and mask low-complexity regions, segmental duplications, and the MHC:

bcftools view -T ^GRCh38_notinalldifficultregions.bed.gz \
  -Oz -o sample.easy.vcf.gz sample.pass.vcf.gz

6. Benchmark the pipeline without stopping at the sample

This is the step most people skip and the one that makes the rest defensible. Run your exact pipeline on a Genome in a Bottle sample with a published truth set (HG002 is the standard, benchmark v4.2.1) and measure precision and recall with hap.py using the vcfeval engine, which handles representational differences in indels correctly:

hap.py HG002_GRCh38_1_22_v4.2.1_benchmark.vcf.gz sample_HG002.norm.vcf.gz \
  -f HG002_GRCh38_1_22_v4.2.1_benchmark_noinconsistent.bed \
  -r GRCh38_no_alt.fna \
  -o happy_out \
  --engine=vcfeval --engine-vcfeval-template GRCh38_sdf \
  --stratification GRCh38_stratification.tsv \
  --threads 16

On 30x PCR-free Illumina with DeepVariant you should see SNV F1 above 0.995 and indel F1 in the 0.99 range inside the high-confidence regions, and considerably worse in low-complexity and segmental-duplication strata. Read the stratified table. Do not rely on the summary line. The summary line is a weighted average dominated by easy genome.

For your own sample, where no truth set exists, use aggregate sanity metrics:

bcftools stats -F GRCh38_no_alt.fna sample.pass.vcf.gz > sample.stats
plot-vcfstats -p plots/ sample.stats

Check that genome-wide Ti/Tv sits around 2.0 to 2.1. A value near 1.5 means false-positive SNVs are diluting the set. Check total variant count against the expected 4.5 to 5 million for a single human genome against GRCh38, and check the het/hom ratio against ancestry expectations (roughly 1.5 to 2.0 for European-ancestry samples, lower for samples with more African ancestry because more of the genome differs from the reference). A sample with unexpectedly high homozygosity may reflect consanguinity or a contaminated library. Run VerifyBamID2 to estimate contamination directly and treat anything above 2 percent as a reason to stop and investigate.

7. Structural variants, copy number, and repeats

Small-variant callers see nothing larger than about 50 bp. Three additional calls are worth making:

# Structural variants (deletions, duplications, inversions, translocations)
configManta.py --bam sample.cram --referenceFasta GRCh38_no_alt.fna \
  --runDir manta_run
manta_run/runWorkflow.py -j 16

# Repeat expansions at known pathogenic loci
ExpansionHunter --reads sample.cram --reference GRCh38_no_alt.fna \
  --variant-catalog variant_catalog_grch38.json \
  --output-prefix sample_eh

Short reads are weak here. Manta’s precision for deletions above 1 kb is good, its precision for insertions is poor, and its ability to resolve breakpoints inside repeats is limited by read length. Expect a few thousand SV calls per genome, most of them common polymorphisms. Long-read sequencing resolves a large class of structural variation that short reads systematically miss, and the gap between what a short-read pipeline calls and what is present in a genome is concentrated in repetitive and structurally complex regions ⁴.

8. Annotate and query

VEP with a local cache is the workhorse. Add gnomAD allele frequencies, ClinVar, and a splicing predictor:

vep -i sample.pass.vcf.gz -o sample.vep.vcf.gz \
  --vcf --compress_output bgzip --cache --offline \
  --assembly GRCh38 --fasta GRCh38_no_alt.fna \
  --everything --pick_allele_gene \
  --custom gnomad.genomes.v4.1.sites.vcf.gz,gnomADg,vcf,exact,0,AF,AF_grpmax,nhomalt \
  --custom clinvar_20250601.vcf.gz,ClinVar,vcf,exact,0,CLNSIG,CLNREVSTAT \
  --plugin SpliceAI,snv=spliceai_scores.snv.hg38.vcf.gz,indel=spliceai_scores.indel.hg38.vcf.gz \
  --fork 16

--pick_allele_gene collapses transcript annotations to one per gene per allele, which keeps the output manageable. If you are looking at a specific gene, drop it and read every transcript, because consequence changes with transcript choice.

Then filter to something small enough to read. slivar expressions are the cleanest way:

slivar expr --vcf sample.vep.vcf.gz --pass-only \
  --info 'INFO.gnomADg_AF < 0.001 && variant.FILTER == "PASS"' \
  --sample-expr 'het:sample.GQ > 20 && sample.DP > 12 && sample.AB > 0.25' \
  -o rare.vcf.gz

A genome yields roughly 100 to 200 protein-truncating variants and tens of thousands of rare coding variants. Almost all of them mean nothing for the individual. Population-scale sequencing established that carrying loss-of-function variants in genes is the normal state. It is not the exception ⁴.

9. Know where the pipeline stops and a clinician starts

Everything above is measurement. Turning a call into a claim about health is a separate discipline with its own standards, and the boundary is here. Three points that matter:

A single-sample research pipeline is not a clinical result. Clinical sequencing workflows require orthogonal confirmation of reportable variants, documented coverage of the regions reported on, and a defined set of genes analyzed ¹. If a variant in your file looks consequential, the next step is a CLIA/CAP laboratory that will re-test it from a fresh sample. A deeper dive into your own BAM is not the next step.

False positives concentrate in predictable places: homopolymer runs longer than 8 bp, short tandem repeats, segmental duplications, and pseudogene-bearing regions. Open the alignment in IGV before believing any single call. Look for reads all in one orientation, all with the same soft-clip position, or a variant supported only by read ends.

Pathogenicity classification is a formal process. It is not a database lookup. A ClinVar “pathogenic” label in a gene with a phenotype you do not have, at an allele frequency inconsistent with the disease prevalence, often reflects a stale submission. Interpretation requires a genetic counselor or medical geneticist.

Common problems

CRAM fails to decode. The reference used at compression time is not the one you have. Check the @SQ M5 checksums against your FASTA with samtools dict, and find the exact matching reference.

Variant counts far above 5 million. Usually contamination, a filtering step that did not apply, or calls in ALT/decoy contigs that should have been excluded. Run bcftools stats per chromosome and see where the excess lives.

Ti/Tv near 1.5 genome-wide. False-positive SNVs. Check that you filtered on PASS, check contamination, and check that you did not accidentally call on a reference with ALT contigs while using a caller configured for a no-ALT set.

Indels disagree between two VCFs of the same sample. The cause is representation. It is not biology. Run bcftools norm -m -any -f ref on both, then compare with rtg vcfeval rather than bcftools isec, because vcfeval compares haplotypes rather than rows.

Chromosome X and Y genotypes look wrong. Most callers assume diploid everywhere unless told otherwise. Set the appropriate ploidy for non-PAR regions, and confirm reported sex against the data before you interpret anything sex-linked.

Results change between runs. Non-deterministic threading. Set -K 100000000 in bwa-mem2, pin the container tag rather than using :latest, and record the reference FASTA checksum with the output. Reproducibility across pipeline implementations is achievable, but only if the steps and parameters are specified rather than assumed ².

Footnotes

1. Daniel C. Koboldt. Best practices for variant calling in clinical sequencing. Genome Medicine, 2020. https://doi.org/10.1186/s13073-020-00791-w ↩ ↩²

2. Allison Regier, Yossi Farjoun, David E. Larson, et al. Functional equivalence of genome sequencing analysis pipelines enables harmonized variant calling across human genetics projects. Nature Communications, 2018. https://doi.org/10.1038/s41467-018-06159-4 ↩ ↩² ↩³

3. Donald Freed, Rafael Aldana, Jessica A. Weber, et al. The Sentieon Genomics Tools - A fast and accurate solution to variant calling from next-generation sequence data. bioRxiv (Cold Spring Harbor Laboratory), 2017. https://doi.org/10.1101/115717 ↩

4. Tuuli Lappalainen, Alexandra J. Scott, Margot Brandt, et al. Genomic Analysis in the Age of Human Genome Sequencing. Cell, 2019. https://doi.org/10.1016/j.cell.2019.02.032 ↩ ↩²