How to Interpret Your Own Variants After the VCF Is Written
Variant interpretation software is the layer that sits between a VCF file and a sentence in English about a specific position in your genome. For a personal whole genome, we would build it as four separate stages rather than buy one platform: annotate with Ensembl VEP running offline against a local cache, filter with bcftools and slivar using population frequency and inheritance logic, score the survivors with CADD and gene-level constraint metrics, then look up each remaining candidate by hand in ClinVar and a guideline-aware aggregator such as VarSome or Franklin. The commercial platforms that rank for this search (Franklin, VarSome, Golden Helix) are useful at the last stage and mostly redundant at the first three. Nothing in this stack produces a diagnosis, and a variant that survives to the end of it is a question for a clinical geneticist, not an answer.
What the stack has to do
A 30x whole genome aligned to GRCh38 produces roughly four to five million positions where you differ from the reference. Almost all of them are common, non-coding, and uninteresting. Perhaps ten to twelve thousand are nonsynonymous coding changes, a few hundred are predicted loss-of-function, and a couple of dozen are rare loss-of-function variants in genes where loss of function is rarely observed in the population. The job of interpretation software is to run that funnel without silently dropping the one variant you cared about, and to attach enough evidence to each survivor that a human can reason about it.
This is a different problem from variant calling, and the two are often confused. Calling asks whether a difference from the reference is real given the reads. Interpretation asks whether a real difference matters. The statistical framing of the second question — what you expect to see by chance in a single genome, and how much power you have to say anything about a variant you have seen once — was laid out early in the exome era and has not changed 1.
Calling, briefly, because interpretation inherits its errors
If you are starting from FASTQ or CRAM, the two defensible choices are GATK HaplotypeCaller in its standard best-practices form and DeepVariant, Google’s convolutional model that calls from pileup images. We prefer DeepVariant for single-sample germline work: it needs no variant quality score recalibration step, it is less sensitive to library idiosyncrasies, and its indel calls in low-complexity regions are noticeably cleaner. Run it as run_deepvariant --model_type=WGS --ref=GRCh38.fa --reads=sample.cram --output_vcf=sample.vcf.gz, budget a few hours on a 32-core machine, and keep the gVCF so you can distinguish “reference” from “no coverage” later. SpeedSeq showed a decade ago that the whole path from reads to annotated variants can run in hours rather than days on commodity hardware, and every modern pipeline is a descendant of that argument 2.
Whichever caller you use, two failure modes propagate straight into interpretation. Regions with no coverage produce absence, not reference, and a pipeline that treats them as reference will quietly tell you that you do not carry a variant you have never been tested for. Segmental duplications and paralogous genes (PMS2, SMN1, HBA1/HBA2, CYP21A2) produce mismapping that yields confident, wrong calls. Short-read data is genuinely unreliable in those regions, and large structural variants and repeat expansions need either long reads or a different physical method entirely 3.
Annotation: VEP, with ANNOVAR as the fast alternative
Annotation attaches biological meaning to coordinates: which transcript, which consequence, which protein position, what the allele frequency is in reference populations. We use Ensembl VEP offline:
vep --offline --cache --dir_cache ~/.vep --assembly GRCh38 \
--fasta GRCh38.fa --everything --pick_allele_gene \
--plugin CADD,whole_genome_SNVs.tsv.gz \
--custom gnomad.genomes.v4.sites.vcf.gz,gnomAD,vcf,exact,0,AF,AF_popmax \
-i sample.vcf.gz -o sample.vep.vcf.gz --vcf --compress_output bgzip
Two flags matter more than the rest. --pick_allele_gene controls which of a gene’s transcripts defines the reported consequence, and the default heuristic will sometimes pick a transcript that is not the one the clinical literature uses. For any variant you plan to act on, check the consequence against the MANE Select transcript explicitly. The --custom line is how you get population frequency in, and gnomAD frequency is the single most powerful filter you have: a variant seen at one percent in any ancestry group is not causing a rare Mendelian disease.
ANNOVAR remains the fastest table-based annotator and is the right choice if you are annotating many samples or want simple tab-separated output rather than an INFO-field bundle 4. The tradeoff is transcript model freshness and the need to maintain your own database directory. VEP’s cache is a single versioned download and stays consistent with Ensembl releases, which is why we default to it.
Filtering: where the funnel is built
After annotation, filtering is a few lines of expression logic. slivar is the tool we reach for because its expressions are JavaScript over the annotated record and are therefore auditable:
slivar expr --vcf sample.vep.vcf.gz --pass-only \
--info 'INFO.gnomAD_AF < 0.001 && variant.FILTER == "PASS"' \
--sample-expr 'het:sample.het && sample.GQ > 20 && sample.DP > 10' \
-o candidates.vcf.gz
Depth and genotype quality thresholds do more work than people expect. Most spurious candidates in a single genome are low-depth heterozygous calls at the edge of a repeat. For a solo genome without parents, you filter on frequency, predicted consequence, and gene-level plausibility. With a parent-child trio, you can add de novo status, which is a far stronger filter, and you can compare the observed de novo count against a mutational model rather than assuming any de novo variant is meaningful 5.
Scoring and knowledge lookup
CADD gives a single integrated deleteriousness score across the whole genome, including non-coding positions, by training on the difference between fixed human-derived alleles and simulated variants 6. Read it on the PHRED-scaled axis: 20 means top one percent of all possible substitutions, 30 means top 0.1 percent. We use it as a ranking heuristic and never as evidence on its own. The same applies to REVEL, AlphaMissense, and SpliceAI, each of which is better than CADD for its specific question and none of which is a classification.
Classification itself follows the ACMG/AMP framework, which combines population frequency, computational prediction, functional data, segregation, and case reports into a five-tier call from benign to pathogenic. Databases now exist that precompute most of the machine-evaluable ACMG criteria genome-wide, which is what makes a lookup tool useful rather than decorative 7. In practice we check every surviving candidate in ClinVar first, read the actual submitter assertions and their review status rather than the aggregate star rating, then cross-check in VarSome or Franklin to see which criteria fired and why. Where those two tools disagree, the disagreement is the finding, and it usually traces to a transcript choice or a stale frequency estimate.
Where this breaks
The dominant result of a personal genome is the variant of uncertain significance: rare, plausibly damaging, in a gene with a known disease association, and supported by no case evidence at all. There is no software that resolves this, and the output is “unknown”. The second failure is ancestry bias in reference databases, which inflates apparent rarity for anyone not of European descent and therefore inflates candidate lists. The third is that whole-genome sequencing generates findings unrelated to whatever you sequenced for, including carrier status and risk alleles that touch relatives who never consented to anything, which is a known ethical problem rather than a technical one 8. Clinical laboratory practice around returning such results is regulated and varies by country 9.
Questions people also ask
Which tool is commonly used for variant calling? GATK HaplotypeCaller is the most widely used and the reference standard for most published pipelines. DeepVariant is the one we would run on a single germline genome, and bcftools mpileup | bcftools call remains fine for quick checks on targeted regions.
How is variant calling done? Reads are aligned to a reference (bwa-mem2 or minimap2), duplicates are marked, and a caller then reassembles reads in each active region and computes genotype likelihoods at positions where the pileup deviates from the reference. The output is a VCF with per-sample genotypes, depth, and quality fields.
What tools do bioinformaticians use for genomics? For short-read germline work: bwa-mem2 or minimap2 for alignment, samtools and bcftools for manipulation, DeepVariant or GATK for calling, Manta or Delly for structural variants, VEP or ANNOVAR for annotation 4, slivar or bcftools filter expressions for filtering, and IGV for looking at the reads when a call seems too good to be true.
How long does sequencing take? Sequencing chemistry on a modern instrument runs one to two days, alignment and calling for a 30x genome take a few hours on a well-provisioned machine, and annotation and filtering take minutes once the caches are local. Interpretation of an ambiguous variant is the slow step and is measured in weeks of literature reading and clinical consultation.
Is sequencing worth it? If you want a permanent, reanalyzable molecular baseline and you are comfortable with most findings being uninterpretable today, yes, because reanalysis against a better database in three years costs nothing beyond compute. If you want a specific clinical answer about a specific condition, a targeted, accredited clinical test ordered through a physician is the faster and more defensible route.
Oak builds longitudinal molecular profiles of individuals: whole-genome sequencing, RNA sequencing, proteomics, blood biomarkers, and continuous glucose data, integrated into one model of you. Build your profile.
Footnotes
-
David B. Goldstein, Andrew Allen, Jonathan Keebler, et al. Sequencing studies in human genetics: design and interpretation. Nature Reviews Genetics, 2013. https://doi.org/10.1038/nrg3455 ↩
-
Colby Chiang, Ryan M Layer, Gregory G Faust, et al. SpeedSeq: ultra-fast personal genome analysis and interpretation. Nature Methods, 2015. https://doi.org/10.1038/nmeth.3505 ↩
-
Michal Levy-Sakin, Yuval Ebenstein. Beyond sequencing: optical mapping of DNA in the age of nanotechnology and nanoscopy. Current Opinion in Biotechnology, 2013. https://doi.org/10.1016/j.copbio.2013.01.009 ↩
-
K. Wang, M. Li, H. Hakonarson. ANNOVAR: functional annotation of genetic variants from high-throughput sequencing data. Nucleic Acids Research, 2010. https://doi.org/10.1093/nar/gkq603 ↩ ↩2
-
James S. Ware, Kaitlin E. Samocha, Jason Homsy, et al. Interpreting de novo Variation in Human Disease Using denovolyzeR. Current Protocols in Human Genetics, 2015. https://doi.org/10.1002/0471142905.hg0725s87 ↩
-
Philipp Rentzsch, Daniela Witten, Gregory M Cooper, et al. CADD: predicting the deleteriousness of variants throughout the human genome. Nucleic Acids Research, 2018. https://doi.org/10.1093/nar/gky1016 ↩
-
Zheng Wang, Guihu Zhao, Zhaopo Zhu, et al. VarCards2: an integrated genetic and clinical database for ACMG-AMP variant-interpretation guidelines in the human whole genome. Nucleic Acids Research, 2023. https://doi.org/10.1093/nar/gkad1061 ↩
-
Holly K. Tabor, Benjamin E. Berkman, Sara Chandros Hull, et al. Genomics really gets personal: How exome and whole genome sequencing challenge the ethical framework of human genetics research. American Journal of Medical Genetics Part A, 2011. https://doi.org/10.1002/ajmg.a.34357 ↩
-
Ros Hastings, Guido de Wert, Brian Fowler, et al. The changing landscape of genetic testing and its impact on clinical and laboratory services and research in Europe. European Journal of Human Genetics, 2012. https://doi.org/10.1038/ejhg.2012.56 ↩