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

By the end of this you will have a bgzipped, tabix-indexed VCF of your small variants (SNVs and indels) called against GRCh38. It will be filtered to PASS with genotype quality thresholds, and annotated with gene consequence and population allele frequency. You will also have a precision/recall number for your own pipeline measured against a truth set. You need: paired-end FASTQ or CRAM from a 30x-or-deeper short-read run (roughly 100 GB of gzipped FASTQ), a machine with 32 cores and 64 GB of RAM, about 1 TB of scratch disk, and Docker or Singularity. Everything here runs on a single box in under a day. Cloud is optional and mostly a cost decision.

1. Get the vocabulary straight

A variant is a position where your sequence differs from the reference genome. A mutation is a variant that arose as a change from a parent sequence: inherited from your parents (germline, present in every cell) or acquired during your life in some cell lineage (somatic, present in a subset of cells). People use the words interchangeably, which is why “mutation calling” and “variant calling” return the same search results. The distinction that matters operationally is germline vs somatic, because they require different callers, different depths, and different filters.

A variant call is a single row in a VCF with these fields.

• chromosome
• position
• reference allele
• alternate allele
• a quality score
• per-sample fields including genotype (0/1 heterozygous, 1/1 homozygous alternate), read depth, and allelic depth

Calling is a statistical inference from noisy alignments. A call is not a readout. Every call has a probability of being wrong, and the whole point of step 7 is to measure yours.

2. Verify the raw data before you spend compute on it

md5sum -c checksums.md5
samtools quickcheck -v *.cram && echo "cram ok"
fastqc -t 8 -o qc/ sample_R1.fastq.gz sample_R2.fastq.gz
seqkit stats sample_R*.fastq.gz

Check three things. Total bases: 30x of a 3.1 Gb genome is ~93 Gbp, so two FASTQs totaling less than ~90 Gbp of sequence are not a 30x run no matter what the invoice says. Per-base quality: a healthy NovaSeq run holds Q30+ across most of the read, with a tail-end drop in the last 10-15 bases. Duplication rate: above ~20% usually means a low-complexity library, and you will lose effective depth after deduplication.

If you received a CRAM instead of FASTQ, confirm which reference it was compressed against. CRAM stores only differences from the reference, so the wrong FASTA gives you garbage or a hard error:

samtools view -H sample.cram | grep '^@SQ' | head -3
samtools view -H sample.cram | grep -o 'UR:[^\t]*' | head -1

3. Choose a reference and stick to it

Use the GRCh38 no-alt analysis set (GCA_000001405.15_GRCh38_no_alt_analysis_set.fna). It includes decoy contigs and has the alt haplotypes excluded, which prevents reads from being split across alt scaffolds and dropping to MAPQ 0. Do not use a plain UCSC hg38 FASTA and do not mix chr1 and 1 naming across your reference, dbSNP, gnomAD, and ClinVar files. Naming mismatches are the single most common reason an annotation step silently produces zero annotations.

T2T-CHM13 is a better genome in every biological sense, and it resolves regions GRCh38 simply does not contain. The problem is that most annotation resources, clinical variant databases, and population frequency panels are still keyed to GRCh38. Our recommendation: call on GRCh38 as your primary, and if you care about specific hard regions, run a second alignment against CHM13 rather than switching wholesale.

4. Align, mark duplicates, and write CRAM

bwa-mem2 index -p ref/GRCh38 ref/GRCh38_no_alt.fna

bwa-mem2 mem -t 32 -K 100000000 -Y \
  -R '@RG\tID:S1_L1\tSM:S1\tLB:lib1\tPL:ILLUMINA\tPU:HXXXXDSX.1.ACGTACGT' \
  ref/GRCh38 sample_R1.fastq.gz sample_R2.fastq.gz \
| samtools fixmate -m -u -@ 4 - - \
| samtools sort -u -@ 8 -m 3G -T /scratch/srt - \
| samtools markdup -@ 8 --reference ref/GRCh38_no_alt.fna \
    -O cram,level=7,seqs_per_slice=10000 - S1.cram

samtools index S1.cram
samtools stats -@ 8 --reference ref/GRCh38_no_alt.fna S1.cram > S1.stats

Two flags earn their keep. -K 100000000 fixes the number of bases processed per batch, which makes output byte-identical regardless of thread count. Without it you cannot reproduce a run exactly, which makes debugging a downstream difference miserable. -Y uses soft clipping for supplementary alignments, which structural-variant callers need.

Piping through fixmate | sort | markdup avoids writing three intermediate 80 GB files. On 32 cores this takes roughly 6-10 hours for 30x. Check samtools stats for the insert size distribution (a single mode around 300-450 bp) and the error rate (typically 0.2-0.5%).

5. Call germline small variants

Use DeepVariant. It is a convolutional classifier over pileup images rather than a hand-built likelihood model, and on standard Illumina WGS it beats GATK HaplotypeCaller on indels by a wide enough margin that the choice is not close for a single-sample germline run.

singularity run -B /data:/data docker://google/deepvariant:1.6.1 \
  /opt/deepvariant/bin/run_deepvariant \
  --model_type=WGS \
  --ref=/data/ref/GRCh38_no_alt.fna \
  --reads=/data/S1.cram \
  --output_vcf=/data/S1.vcf.gz \
  --output_gvcf=/data/S1.g.vcf.gz \
  --num_shards=32 \
  --intermediate_results_dir=/scratch/dv_tmp

Expect 3-6 hours on 32 cores. Always emit the gVCF as well as the VCF. The gVCF records reference confidence in non-variant blocks, which is what lets you distinguish “no variant here” from “no data here” later. That distinction matters enormously when you go looking for a specific position and find nothing.

Keep GATK HaplotypeCaller in reserve for two cases. The first is joint calling across a family or cohort (GenotypeGVCFs still has real advantages there). The second is any workflow where you need a documented, deterministic, non-neural-net model. Large cohort projects lean on exactly this kind of configuration-driven, reproducible pipeline structure, and that reproducibility is worth more than a few tenths of a percent of F1 when you are comparing samples to each other.¹

Mitochondrial variants need a separate path. Numts (nuclear copies of mtDNA) confuse standard callers, and mtDNA heteroplasmy is a continuous allele fraction. Heteroplasmy is not a genotype. Run GATK’s mitochondria mode against the chrM-only reference with --mitochondria-mode.

6. Open and read the VCF

A VCF is a tab-delimited text file with ## metadata headers, one #CHROM column header line, and one row per variant site. It is nearly always bgzip-compressed with a .tbi index, so cat and less on the raw file give you binary. Use bcftools.

bcftools view -h S1.vcf.gz | grep '^##FORMAT'
bcftools view S1.vcf.gz chr7:117,480,000-117,670,000 | less -S

bcftools query \
  -f '%CHROM\t%POS\t%REF\t%ALT\t%QUAL\t%FILTER[\t%GT\t%DP\t%AD\t%GQ]\n' \
  -r chr17:43044295-43125483 S1.vcf.gz | head

Reading a row: GT of 0/1 means one reference and one alternate allele. DP is the read depth at that site after filtering. AD is the per-allele depth, e.g. 18,16 for a clean het. GQ is the phred-scaled confidence in the genotype call, so GQ 20 is a 1% chance the genotype is wrong and GQ 50 is 1 in 100,000. A heterozygous call with AD of 28,3 is a likely false positive, mapping artifact, or low-level somatic event. Such a call is not a germline het.

Sanity-check the whole file before you trust any single row:

bcftools stats -F ref/GRCh38_no_alt.fna S1.vcf.gz > S1.vchk
plot-vcfstats -p plots/ S1.vchk

For a 30x human WGS you should see roughly 4.5-5.0 million SNVs and 800,000-1,000,000 indels, with a transition/transversion ratio near 2.0-2.1 genome-wide. A Ti/Tv of 1.6 means your filters are letting sequencing noise through. Counts far above 5.5 million usually mean contamination from a second individual.

7. Filter, then annotate

bcftools view -f PASS -i 'FMT/GQ>=20 && FMT/DP>=10 && FMT/DP<=100' \
  -Oz -o S1.pass.vcf.gz S1.vcf.gz
bcftools index -t S1.pass.vcf.gz

The upper depth bound removes collapsed repeats and centromeric pileups, where apparent depth of 300x reflects many paralogous regions mapping to one place.

vep -i S1.pass.vcf.gz -o S1.vep.vcf.gz \
  --cache --offline --dir_cache ~/.vep --assembly GRCh38 \
  --vcf --compress_output bgzip --fork 8 --buffer_size 20000 \
  --everything --pick_allele_gene \
  --custom file=gnomad.genomes.v4.1.sites.vcf.bgz,short_name=gnomADg,format=vcf,type=exact,fields=AF_grpmax \
  --plugin CADD,snv=whole_genome_SNVs.tsv.gz,indels=gnomad.genomes.indel.tsv.gz

Then narrow to things worth a second look:

filter_vep -i S1.vep.vcf.gz \
  -filter "gnomADg_AF_grpmax < 0.001 and (IMPACT is HIGH or CADD_PHRED > 25)" \
  -o S1.rare_impactful.vcf

Two notes on allele frequency filtering. First, AF_grpmax (max frequency across ancestry groups) is the right field, because a variant that is 3% in one population and absent from the aggregate is common. Such a variant is not rare. Second, if your ancestry is underrepresented in gnomAD, a genuinely common variant in your population can look novel. Deep sequencing of population-specific cohorts repeatedly finds large numbers of variants absent from existing reference panels, which is precisely why ancestry-matched frequency data changes what survives a rarity filter.²

Nothing produced here is a diagnosis. A research-grade pipeline on a research-grade sample produces hypotheses. Any variant you intend to act on needs orthogonal confirmation in a CLIA/CAP-certified laboratory and interpretation with a clinical geneticist or genetic counselor, who will weigh it against your family history and phenotype in ways no annotation pipeline does.

8. Benchmark your pipeline against a truth set

This step is what separates a pipeline you can reason about from one you are guessing at. Run your exact commands on public HG002 FASTQs (Genome in a Bottle, Ashkenazi son), then compare to the truth VCF.

hap.py HG002_truth.vcf.gz S1_HG002.vcf.gz \
  -f HG002_confident_regions.bed \
  -r ref/GRCh38_no_alt.fna \
  --engine=vcfeval --engine-vcfeval-template ref/GRCh38.sdf \
  --threads 16 -o bench/hg002

For 30x Illumina plus DeepVariant inside the GIAB high-confidence regions, SNV F1 above 0.995 and indel F1 above 0.99 are normal. If you are well below that, something is wrong with your alignment or filtering. The caller is not the problem.

The more interesting question is what happens outside those regions. The high-confidence BED deliberately excludes segmental duplications, long tandem repeats, and other regions where the truth itself was uncertain. Benchmarks built from complete diploid assemblies now extend truth into much of that excluded territory, and accuracy in the newly covered regions is substantially lower than the headline numbers suggest.³ Run your benchmark twice, once inside the classic confident regions and once against the expanded set. Treat the second number as the real one when you go looking at a gene that sits in a duplicated region (PMS2, SMN1, CYP21A2 or HBA1/HBA2).

9. If you are calling somatic mutations

Germline calling asks: which of two alleles did you inherit. Somatic calling asks: is this alternate allele present in some fraction of cells and absent from the germline. That is a different problem, and it needs a matched normal.

gatk Mutect2 -R ref/GRCh38_no_alt.fna \
  -I tumor.cram -I normal.cram -normal NORMAL_SM \
  --germline-resource af-only-gnomad.hg38.vcf.gz \
  --panel-of-normals pon.vcf.gz \
  --f1r2-tar-gz f1r2.tar.gz \
  -O unfiltered.vcf.gz

gatk LearnReadOrientationModel -I f1r2.tar.gz -O read-orientation-model.tar.gz
gatk GetPileupSummaries -I tumor.cram -V small_exac_common_3.hg38.vcf.gz \
  -L small_exac_common_3.hg38.vcf.gz -O pileups.table
gatk CalculateContamination -I pileups.table -O contamination.table
gatk FilterMutectCalls -R ref/GRCh38_no_alt.fna -V unfiltered.vcf.gz \
  --contamination-table contamination.table \
  --ob-priors read-orientation-model.tar.gz \
  -O somatic.filtered.vcf.gz

Three things to internalize. Depth determines what you can see: at 30x tumor coverage you are reliably detecting variant allele fractions above roughly 10-15%, which means clonal events only. Subclonal populations need 100-300x, and circulating tumor DNA needs duplex UMI chemistry on top of that. Sample purity and heterogeneity compound the problem, and callers vary considerably in how they handle mixtures of clones at different fractions.⁴ Caller choice matters more than it does for germline: head-to-head comparisons of somatic callers on the same data show meaningfully different sensitivity and specificity profiles, with concordance between callers falling off sharply at low allele fractions.⁵ Running two callers and taking the intersection is a defensible default, at the cost of sensitivity.

Serial sampling is where somatic WGS gets genuinely informative. Sequencing the same person’s tumor at multiple timepoints tracks which subclones expand and which contract, which no single snapshot can show.⁶ Cancer whole-genome sequencing also picks up structural rearrangements and mutational signatures that exome sequencing misses entirely, because most of the informative passenger mutations are not in coding sequence.⁷

One interpretive caution: somatic mutation density is not uniform across the genome. Mismatch repair activity varies regionally, so late-replicating heterochromatin accumulates mutations at several times the rate of early-replicating, actively transcribed regions.⁸ A cluster of somatic calls in a gene means much less if that gene sits in a high-baseline-rate region.

Common problems

The VCF has 8 million variants. Sample contamination from a second individual. Run verifyBamID2 or somalier and check the FREEMIX estimate. Above 0.02 you have a problem.

Zero annotations after VEP. Contig naming mismatch (chr1 vs 1) between your VCF and the cache, or an assembly mismatch (GRCh37 cache against a GRCh38 VCF). Check with bcftools view -h | grep contig.

A gene you care about has no calls. Check coverage before concluding absence: samtools depth -a -r chr7:117479963-117668665 S1.cram | awk '$3<10'. GC-extreme first exons and segmental duplications routinely drop below callable depth. The gVCF tells you which case you are in.

Indel calls pile up in homopolymers. Runs of 8+ identical bases produce polymerase slippage in the library prep, and short reads cannot resolve them. Filter indels in homopolymer tracts of length ≥8 unless you have long reads.

CRAM reads fail or produce nonsense. The reference used for compression is not the reference you are decompressing with. Pass it explicitly: samtools view -T ref.fna.

Two runs of the same pipeline give different VCFs. Missing -K in bwa-mem2, or a thread-count-dependent step downstream. Pin your container digests. Tags alone are not enough.

Tumor-in-normal contamination. If the matched normal is a blood draw from someone with a hematologic malignancy or high circulating tumor DNA, somatic calls get filtered away as germline. Use a tissue normal where possible, and inspect the normal_artifact filter counts.

Footnotes

1. Manuel Rueda, Dietmar Fernandez Orth, Ivo G Gut. CBIcall: a configuration-driven framework for variant calling in large sequencing cohorts. Bioinformatics Advances, 2026. https://doi.org/10.1093/bioadv/vbag232 ↩

2. Masao Nagasaki, Jun Yasuda, Fumiki Katsuoka, et al. Rare variant discovery by deep whole-genome sequencing of 1,070 Japanese individuals. Nature Communications, 2015. https://doi.org/10.1038/ncomms9018 ↩

3. Nancy F. Hansen, Nathan Dwarshuis, Hyun Joo Ji, et al. A complete diploid human genome benchmark for personalized genomics. Cell, 2026. https://doi.org/10.1016/j.cell.2026.06.016 ↩

4. Shuangxi Ji, Tong Zhu, Ankit Sethia, et al. Accelerated somatic mutation calling for whole-genome and whole-exome sequencing data from heterogenous tumor samples. Genome Research, 2024. https://doi.org/10.1101/gr.278456.123 ↩

5. Huilei Xu, John DiCarlo, Ravi Vijaya Satya, et al. Comparison of somatic mutation calling methods in amplicon and whole exome sequence data. BMC Genomics, 2014. https://doi.org/10.1186/1471-2164-15-244 ↩

6. Anna Schuh, Jennifer Becq, Sean Humphray, et al. Monitoring chronic lymphocytic leukemia progression by whole genome sequencing reveals heterogeneous clonal evolution patterns. Blood, 2012. https://doi.org/10.1182/blood-2012-05-433540 ↩

7. Hidewaki Nakagawa, Christopher P. Wardell, Mayuko Furuta, et al. Cancer whole-genome sequencing: present and future. Oncogene, 2015. https://doi.org/10.1038/onc.2015.90 ↩

8. Fran Supek, Ben Lehner. Differential DNA mismatch repair underlies mutation rate variation across the human genome. Nature, 2015. https://doi.org/10.1038/nature14173 ↩