How to interpret your own whole-genome variants

By the end of this you will have four artifacts.

• a coverage report that tells you which genes your sequencing can and cannot speak to
• a normalized and annotated VCF with gnomAD frequencies, ClinVar assertions, and in-silico predictions attached to every alternate allele
• a tiered shortlist of maybe 20-80 variants worth human attention
• a pharmacogenomic diplotype table

You need your CRAM or BAM plus the germline VCF, the exact reference FASTA the data was aligned to, roughly 300 GB of free disk (the VEP cache and gnomAD sites files dominate), 16-32 GB RAM, and Docker. Everything below is free software and public data. Budget an afternoon for setup and an evening for the first pass.

A calibration note before the commands. A single human genome differs from the reference at roughly 4-5 million sites, which is the scale the first genome-wide SNP maps anticipated when they catalogued 1.42 million SNPs at about one per 1.9 kb ¹. Almost all of those differences are common, shared with millions of people, and uninterpretable in isolation. The job is not to explain your genome. The job is to find the handful of sites where the evidence is strong enough to mean something, and about the rest. Nothing here is a diagnosis. If a variant survives your filters and lands in a gene with clinical management implications, that is the point at which you stop and take it to a clinical geneticist or a genetic counselor for confirmatory testing on an independent sample.

1. Establish what build and what file you have

Half of self-interpretation failures start here. Check the reference:

samtools view -H sample.cram | grep -m5 '^@SQ'
bcftools view -h sample.vcf.gz | grep -E '##reference|##contig=<ID=chr1,'

If contigs are named chr1 you are on a UCSC-style GRCh38 or hg19 build. If they are 1, likely Ensembl/1000G GRCh37 or GRCh38. GRCh38 comes in flavors: with or without ALT contigs, with or without the HLA decoys, and with or without the masked false duplications in the analysis set. Mixing a no-ALT reference with an ALT-aware BAM changes mapping quality in exactly the segmental duplications you care about. Record the MD5 of your FASTA and reuse that file for every downstream step.

Then sanity-check the variant set:

bcftools stats sample.vcf.gz | grep -E '^SN|^TSTV'

For a 30x PCR-free WGS you expect roughly 4.0-5.0 million SNVs, 500,000-900,000 indels, and a transition/transversion ratio near 2.0-2.1. A Ti/Tv of 1.6 means you are looking at noise. A SNV count of 3.2 million means the VCF was filtered or you are on a different build than you think.

2. Measure coverage before you trust any negative result

An absent variant call and a region with zero reads look identical in a VCF. Compute callable fraction over the genes you care about:

mosdepth --threads 4 --by genes.bed --fasta GRCh38.fa \
  --mapq 20 --quantize 0:1:10:20: sample sample.cram

The --quantize bins give you a BED of regions at 0x, 1-9x, 10-19x, and 20x+. Intersect the sub-10x bin with your gene panel. Short reads reliably lose coverage in GC-rich first exons, in homopolymers, and in segmental duplications. Those duplications include SMN1/SMN2, CYP2D6/CYP2D7, PMS2/PMS2CL, HBA1/HBA2, STRC and IKBKG. If one of those is on your list, a standard short-read pipeline is not the right instrument and long reads or a targeted assay are. Write the callable fraction per gene into a table and keep it next to your results. Every negative finding should be read as “no variant called in the 94% of this gene that was callable.”

3. Normalize, decompose, and left-align

Annotation databases are keyed on exact CHROM-POS-REF-ALT. Multi-allelic records and right-aligned indels will silently miss ClinVar matches.

bcftools norm -m -any -f GRCh38.fa -c w -Ou sample.vcf.gz \
  | bcftools norm -d exact -Oz -o sample.norm.vcf.gz
bcftools index -t sample.norm.vcf.gz

The -m -any option splits multi-allelics into one row per alternate allele, and the reference FASTA drives left-alignment. The -c w option warns on REF mismatches instead of failing silently, and -d exact drops duplicates created by the split. Compare record counts before and after. A large jump means your caller emitted many multi-allelic sites, which is normal.

4. Annotate with VEP plus the plugins that matter

Run Ensembl VEP offline against a local cache. The plugin set is the part people skip and then regret.

vep --input_file sample.norm.vcf.gz --format vcf --vcf \
  --cache --offline --dir_cache $VEP_CACHE --assembly GRCh38 \
  --fasta GRCh38.fa --fork 8 --buffer_size 20000 \
  --everything --pick_allele_gene --mane_select \
  --plugin AlphaMissense,file=AlphaMissense_hg38.tsv.gz \
  --plugin SpliceAI,snv=spliceai_scores.raw.snv.hg38.vcf.gz,indel=spliceai_scores.raw.indel.hg38.vcf.gz \
  --plugin CADD,snv=whole_genome_SNVs.tsv.gz,indels=gnomad.genomes.indel.tsv.gz \
  --plugin LoFtool --plugin NMD \
  --custom file=clinvar.vcf.gz,short_name=ClinVar,format=vcf,type=exact,fields=CLNSIG%CLNREVSTAT%CLNDN%CLNHGVS \
  --custom file=gnomad.genomes.v4.1.sites.vcf.bgz,short_name=gnomADg,format=vcf,type=exact,fields=AF%AF_nfe%AF_afr%AF_eas%AF_sas%AF_amr%nhomalt \
  --compress_output bgzip --output_file sample.vep.vcf.gz --stats_file vep.html

Notes on the choices. --pick_allele_gene collapses the transcript explosion to one consequence per allele per gene, and --mane_select biases that pick toward the MANE transcript, which is what ClinVar and most labs report against. Keep --everything so you retain HGVS notation, existing variation IDs, and per-transcript consequence terms. Pull ClinVar as a --custom annotation rather than relying on the cache, because ClinVar changes weekly and you want to know the release date of the copy you used. Record it.

SpliceAI and AlphaMissense do different jobs. SpliceAI delta scores flag splice disruption up to 50 bp from the variant, which catches synonymous and deep-intronic variants that consequence terms label as harmless. AlphaMissense gives a pathogenicity class for missense variants where you would otherwise have nothing but conservation. Both are predictions. Neither is evidence sufficient on its own to call a variant pathogenic, and both are used in the ACMG/AMP framework only as supporting-strength evidence.

5. Flatten to a table and tier it

bcftools +split-vep sample.vep.vcf.gz -d -A tab \
  -f '%CHROM\t%POS\t%REF\t%ALT\t[%GT]\t%SYMBOL\t%HGVSc\t%HGVSp\t%Consequence\t%IMPACT\t%gnomADg_AF\t%gnomADg_nhomalt\t%ClinVar_CLNSIG\t%ClinVar_CLNREVSTAT\t%am_class\t%CADD_PHRED\t%SpliceAI_pred_DS_AG\t%SpliceAI_pred_DS_AL\t%SpliceAI_pred_DS_DG\t%SpliceAI_pred_DS_DL\n' \
  -i 'FILTER="PASS"' > sample.flat.tsv

Now build tiers in whatever you already use (pandas, DuckDB, awk). We use three.

Tier 1, ClinVar-driven: CLNSIG contains Pathogenic or Likely_pathogenic, and CLNREVSTAT is criteria_provided,_multiple_submitters,_no_conflicts or better (two stars and up). Two stars is the floor we accept without reading primary literature. One-star and no-assertion records include a meaningful number of entries that later get reclassified.

Tier 2, frequency plus effect: gnomAD global AF below 0.001 and IMPACT=HIGH, and not flagged by LOFTEE as low-confidence. The high-impact classes are stop-gained, frameshift, canonical splice and start-loss. This tier is dominated by false positives in the last exon, in genes with no disease association, and in repeat-rich sequence. Expect 100-300 rows genome-wide before you apply a gene list.

Tier 3, prediction-driven: AF below 0.0001, missense with AlphaMissense likely_pathogenic, or any variant with a SpliceAI delta score above 0.5 regardless of consequence term. This tier is where you find things that are interesting and unresolvable.

Restrict all three to a gene list before you read them. Clinical WGS pipelines drive filtering from phenotype rather than scanning all genes, because the unfiltered genome produces far more plausible-looking variants than any person can adjudicate ². If you have no phenotype to work from, use the ACMG secondary-findings gene list, which is deliberately limited to genes with actionable management. Add a carrier panel if reproductive planning is the question.

6. Read a variant the way a lab would

For each survivor, gather five things before forming a view.

Allele frequency in your own ancestry. The global number is not the figure to use. A variant at 0.4% in East Asian populations and 0.001% globally is not rare in the relevant sense. Population databases are also incomplete outside the ancestries they sampled: deep whole-genome sequencing of 1,070 Japanese individuals recovered large numbers of variants absent from existing reference panels, which is the general shape of the problem ³. Absence from gnomAD is weak evidence when gnomAD has thin coverage of your ancestry.

Homozygote count. nhomalt above zero in gnomAD for a variant claimed to cause a severe recessive condition is a strong argument against the claim.

ClinVar review status and the submitters. Open the variant page. A single submitter asserting pathogenic in 2013 with no evidence codes is different from an expert panel review.

Zygosity and phase. 0/1 in a recessive gene means nothing without a second hit. If you have two candidate variants in one gene, phase them:

whatshap phase --reference GRCh38.fa -o phased.vcf.gz sample.norm.vcf.gz sample.cram

Read-backed phasing works when the two variants are within a read or fragment length. Beyond that you need parental samples or long reads.

Read support. Load the region in IGV. Look at allele balance (a real het sits near 0.5) and at mapping quality. Look at whether the supporting reads all come from one strand, and at whether the site sits in a homopolymer. A surprising fraction of exciting-looking indels evaporate on inspection.

7. Run pharmacogenomics separately

Star-allele calling is a different problem from single-variant interpretation, because the unit is a haplotype and several genes sit in duplicated sequence. Use PharmCAT rather than hand-reading CYP2C19 variants:

pharmcat_vcf_preprocessor -vcf sample.norm.vcf.gz -refFna GRCh38.fa -o pgx/
pharmcat -vcf pgx/sample.preprocessed.vcf.bgz -reporterJson -o pgx/

PharmCAT emits diplotypes and the matching CPIC guideline text. Read CYP2D6 results with suspicion: the gene has hybrid alleles and copy-number variation that short reads call poorly, and PharmCAT delegates it to an external caller. Do not start, stop, or change any medication based on this output. Pharmacogenomic results change prescribing decisions only through a clinician, and clinical labs re-test on a fresh sample before acting.

8. Structural variants, repeats, and the noncoding genome

Your SNV/indel VCF misses whole classes of variation. Two extra runs are worth the compute:

configManta.py --bam sample.cram --referenceFasta GRCh38.fa --runDir manta/
manta/runWorkflow.py -j 8

ExpansionHunter --reads sample.cram --reference GRCh38.fa \
  --variant-catalog variant_catalog_hg38.json --output-prefix eh/sample

Manta calls deletions, duplications, insertions, inversions, and translocations from paired-end and split-read signal. ExpansionHunter genotypes known repeat loci. Both have high false-positive rates on short reads, and both need IGV confirmation of the breakpoint reads. The variety of rearrangement types that whole-genome sequencing exposes, and how badly exome or panel data represents them, is clearest in cancer genome studies where structural variation dominates ⁴. That work is somatic, and your germline profile says nothing about tumor mutations, but the technical lesson about what short-read SV calling can and cannot see carries over.

The noncoding fraction is the limit of this exercise. Regulatory elements, transcribed non-coding RNAs, and chromatin features cover far more of the genome than protein-coding exons do ⁵, and variant interpretation in those regions has no ClinVar-scale evidence base. SpliceAI is the one prediction we trust enough to act on in noncoding sequence, and only at high delta scores near canonical junctions.

9. Version everything and re-run on a schedule

Your VCF is static. The annotations are not. ClinVar reclassifies variants continuously, gnomAD adds cohorts, and predictors get retrained. Pin the versions in a manifest:

reference: GRCh38_full_analysis_set_plus_decoy_hla.fa  md5 64b32de2fc934679c16e83a2bc072064
vep: 112 / cache 112_GRCh38
clinvar: clinvar_20260901.vcf.gz
gnomad: v4.1 genomes
alphamissense: hg38 2023-09
spliceai: 1.3 raw scores

Re-run step 4 onward every six months against fresh ClinVar and diff the Tier 1 table. Reclassification in both directions is common, and a variant that moved from VUS to likely pathogenic is exactly the thing you want to notice.

Common problems

Reference mismatch between VCF and annotation files. Symptom: near-zero ClinVar matches, or VEP warning about REF allele mismatches. Fix: confirm build and contig naming on every input, and use bcftools annotate --rename-chrs rather than hand-editing.

Forgetting to normalize. Symptom: a known pathogenic indel in your data that ClinVar does not match. Indels have multiple valid representations, and only the left-aligned one hits the database key.

Reading a negative as reassurance. If a gene is 82% callable at 20x, you have not ruled anything out. Keep the coverage table beside the results.

Treating CADD or AlphaMissense as a verdict. These scores separate populations of variants. A score does not settle an individual variant. A CADD PHRED of 28 on a missense variant in a gene with no disease link tells you almost nothing.

One-star ClinVar entries. Older submissions from single labs carry meaningful reclassification rates. Two stars is the working floor, and expert-panel or practice-guideline entries are the only ones we treat as settled.

Polygenic scores computed on one sample. plink2 --score will produce a number for a single VCF, but the number is meaningless without a reference cohort of the same ancestry to convert it into a percentile, and PRS weights transfer poorly across ancestries. If you compute one, compute it against 1000 Genomes samples processed through your exact pipeline, and report the percentile with the ancestry group named.

Assuming a germline VCF says anything about somatic mutation. It does not. Different sample, different pipeline, different interpretation rules.

Acting alone on an actionable finding. A pathogenic variant in a cancer-predisposition or cardiomyopathy gene from a research-grade pipeline needs orthogonal confirmation in a clinical laboratory and a conversation with a genetic counselor before it changes anything. Your pipeline produces a hypothesis with a good prior. That is where the technical work stops.

Footnotes

1. The International SNP Map Working Group, Cold Spring Harbor Laboratories:, Ravi Sachidanandam, et al. A map of human genome sequence variation containing 1.42 million single nucleotide polymorphisms. Nature, 2001. https://doi.org/10.1038/35057149 ↩

2. Frederik Otzen Bagger, Line Borgwardt, Andreas Sand Jespersen, et al. Whole genome sequencing in clinical practice. BMC Medical Genomics, 2024. https://doi.org/10.1186/s12920-024-01795-w ↩

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

4. Australian Pancreatic Cancer Genome Initiative, Nicola Waddell, Marina Pajic, et al. Whole genomes redefine the mutational landscape of pancreatic cancer. Nature, 2015. https://doi.org/10.1038/nature14169 ↩

5. Ewan Birney, Paul Flicek, Damian Keefe, et al. Identification and analysis of functional elements in 1% of the human genome by the ENCODE pilot project. Nature, 2007. https://doi.org/10.1038/nature05874 ↩