How to Search Your Own Variants: From VCF to a Shortlist You Can Read

By the end of this you will have a normalized, indexed VCF of your own genome annotated with ClinVar significance, gnomAD allele frequencies, and VEP consequence terms, plus a filter chain that takes it from roughly 4.5 million variants down to a shortlist of a few dozen you can read one at a time in a browser. You need: your WGS VCF (or gVCF) and the matching reference FASTA, about 100 GB of free disk, 16 GB of RAM, and bcftools 1.19+, tabix, samtools, and either Ensembl VEP with an offline cache or a Docker runtime to host it. Everything below assumes GRCh38. If your data is on hg19 or T2T-CHM13, read step 1 before running anything, because coordinates are the single most common source of wrong answers in variant work ¹.

1. Find out exactly what you have

Before any annotation, establish three facts: the reference build, the contig naming convention, and whether the file is a filtered callset or a gVCF with reference blocks.

bcftools view -h your.vcf.gz | grep -E '^##(reference|contig=<ID=(chr)?1,)'
bcftools view -h your.vcf.gz | grep -c '^##contig'
bcftools stats your.vcf.gz > stats.txt
grep -E '^SN' stats.txt

The SN block gives counts. For a 30x PCR-free short-read genome called against GRCh38, expect roughly 4.3–5.0 million SNPs and 800,000–1,000,000 indels for a single sample, with the total varying by ancestry (non-African genomes sit at the lower end). Titransversion ratio (bcftools stats reports it as TSTV) should be near 2.0–2.1 genome-wide and 2.8–3.0 in coding regions. A Ti/Tv of 1.6 means your callset is loaded with false positives and every downstream filter will be polluted.

If contig=<ID=chr1 appears, you have UCSC-style names. ClinVar’s VCF uses Ensembl-style (1, 2 and so on through MT). You will need a rename map either way:

for i in $(seq 1 22) X Y MT; do echo -e "$i\tchr$i"; done > ens2ucsc.txt
bcftools annotate --rename-chrs ens2ucsc.txt -Oz -o clinvar.chr.vcf.gz clinvar.vcf.gz

If the file has <NON_REF> alleles and END tags, it is a gVCF. Convert first with bcftools view -e 'ALT="<NON_REF>"' after genotyping, or ask whoever produced it for the sample VCF. Also check the caller: DeepVariant, GATK HaplotypeCaller, and graph-based callers disagree meaningfully on indels in low-complexity and repeat-rich regions, and a nontrivial share of true variation sits in exactly those regions ^{2 3}.

2. Normalize, left-align, and deduplicate

This step is non-negotiable. An indel represented in two different but equivalent ways will not match the ClinVar record for the same indel, and you will get a silent false negative on a variant that is in the database.

bcftools norm -m -any -f GRCh38.fa -Ou your.vcf.gz \
  | bcftools norm -d exact -Oz -o norm.vcf.gz -W=tbi

-m -any splits multiallelic sites into one row per ALT. -f with the reference left-aligns and trims indels. -d exact drops duplicate records created by the split. -W=tbi writes the index in the same pass (older bcftools: run tabix -p vcf norm.vcf.gz after).

Verify the reference matches. If your FASTA is a different GRCh38 patch or includes/excludes ALT contigs, bcftools norm prints REF_MISMATCH warnings to stderr. More than a handful means the wrong FASTA:

bcftools norm -c w -f GRCh38.fa norm.vcf.gz 2>&1 >/dev/null | head

If your VCF is hg19 and you want to use current ClinVar and gnomAD, lift over rather than downgrade the databases. CrossMap vcf hg19ToHg38.over.chain GRCh38.fa in.vcf.gz out.vcf, then re-normalize. Liftover drops 0.1–0.5% of sites and silently flips strand on some. The dropped file is written out, so read it. Cross-build workflows are common enough that platforms now handle T2T, hg38, and hg19 side by side ¹.

3. Annotate with ClinVar

Pull the weekly release and attach the fields you care about.

wget https://ftp.ncbi.nlm.nih.gov/pub/clinvar/vcf_GRCh38/clinvar.vcf.gz{,.tbi}

bcftools annotate \
  -a clinvar.chr.vcf.gz \
  -c INFO/CLNSIG,INFO/CLNREVSTAT,INFO/CLNDN,INFO/CLNHGVS,INFO/ALLELEID \
  -Oz -o norm.clinvar.vcf.gz -W=tbi norm.vcf.gz

CLNREVSTAT is the field most people skip and it matters more than CLNSIG. A “Pathogenic” call with no_assertion_criteria_provided is one submitter’s opinion from 2011. A “Pathogenic” call with reviewed_by_expert_panel or practice_guideline carries real weight. Count what you got:

bcftools query -f '%INFO/CLNSIG\t%INFO/CLNREVSTAT\n' norm.clinvar.vcf.gz \
  | grep -v '^\.' | sort | uniq -c | sort -rn | head -30

A typical genome matches 10,000–30,000 ClinVar records. The overwhelming majority are Benign, Likely_benign, or Uncertain_significance. Expect somewhere between 2 and 10 hits labeled Pathogenic or Likely_pathogenic with at least one-star review status, and most of those will be recessive carrier alleles for conditions you do not have.

4. Add population frequencies

Allele frequency is the strongest single filter you have. Anything present in 1% of gnomAD is not a rare high-penetrance disease allele, whatever ClinVar says about it.

gnomAD v4 exomes+genomes is ~2.4 TB of VCF, which you do not want. Use echtvar, which encodes selected INFO fields into a compressed zip archive (~2–4 GB for the AF fields) and annotates a whole genome in under a minute:

echtvar encode gnomad_v4.zip conf.json gnomad.genomes.v4.*.sites.vcf.bgz
echtvar anno -e gnomad_v4.zip norm.clinvar.vcf.gz norm.clinvar.gnomad.bcf

The conf.json names the fields to carry over (AF, AF_popmax/AF_grpmax, nhomalt, and AC). If you would rather not build the archive, pre-built gnomAD echtvar files and slivar gnotate files are published and are the faster path.

Add nhomalt. A variant with dozens of homozygotes in gnomAD is very unlikely to be a severe recessive allele, and that single number resolves a lot of ambiguous ClinVar entries on its own.

5. Get functional consequences from VEP

ClinVar only knows about variants someone has submitted. For the other 4.4 million, you need predicted consequence.

vep -i norm.clinvar.gnomad.bcf \
  -o vep.vcf.gz --vcf --compress_output bgzip \
  --cache --dir_cache $HOME/.vep --offline --assembly GRCh38 \
  --fasta GRCh38.fa \
  --everything --pick_allele_gene --canonical --mane \
  --fork 8 --buffer_size 50000 \
  --plugin CADD,whole_genome_SNVs.tsv.gz,gnomad.genomes.indel.tsv.gz \
  --plugin SpliceAI,snv=spliceai_scores.raw.snv.hg38.vcf.gz,indel=spliceai_scores.raw.indel.hg38.vcf.gz

On 8 cores this takes 40–90 minutes for a whole genome. Notes on the flags:

• --pick_allele_gene emits one consequence per allele per gene instead of one per transcript. Without it, a single variant can produce 40 CSQ entries and your downstream parsing gets ugly. Use --mane so the chosen transcript is the MANE Select one where it exists, which is what clinical labs and HGVS reports use.
• --everything turns on SIFT, PolyPhen, regulatory features, and the frequency fields. It roughly doubles runtime.
• SpliceAI is the plugin worth the disk space. A synonymous or deep-intronic variant with a delta score above 0.5 is a real candidate that consequence terms alone would have thrown away.

Then flatten to a table with +split-vep:

bcftools +split-vep vep.vcf.gz -d -f \
 '%CHROM\t%POS\t%REF\t%ALT\t%SYMBOL\t%Consequence\t%HGVSc\t%HGVSp\t%CADD_PHRED\t%SpliceAI_pred_DS_AG\t%INFO/gnomad_AF_grpmax\t%INFO/CLNSIG\t%INFO/CLNREVSTAT\t[%GT]\n' \
 -i 'GT!="0/0"' > variants.tsv

Load variants.tsv into DuckDB or pandas and the rest is SQL.

6. Build the filter chain

Two separate queries, because they answer different questions.

Query A, the known-variant question: what do curated databases already say about alleles I carry?

SELECT * FROM variants
WHERE CLNSIG LIKE '%athogenic%'
  AND CLNSIG NOT LIKE '%Conflicting%'
  AND CLNREVSTAT NOT LIKE '%no_assertion%'
  AND (gnomad_AF_grpmax < 0.01 OR gnomad_AF_grpmax IS NULL)
ORDER BY CLNREVSTAT DESC;

Expect 2–15 rows. Read every one.

Query B, the unknown-variant question: which variants look damaging even though nobody has curated them?

SELECT * FROM variants
WHERE (gnomad_AF_grpmax < 0.001 OR gnomad_AF_grpmax IS NULL)
  AND nhomalt < 5
  AND (
    Consequence LIKE '%stop_gained%' OR
    Consequence LIKE '%frameshift%' OR
    Consequence LIKE '%splice_acceptor%' OR
    Consequence LIKE '%splice_donor%' OR
    (Consequence LIKE '%missense%' AND CADD_PHRED > 25) OR
    SpliceAI_DS_max > 0.5
  );

That lands around 100–400 rows genome-wide. To get to something readable, intersect with a gene list rather than loosening the thresholds. Options are the ACMG SF v3.2 secondary-findings gene list (81 genes), a PanelApp green-gene panel for a phenotype you have reason to care about, or the PharmGKB very-important-pharmacogene set ⁴. Gene-list intersection is how family-based WGS prioritization studies get from millions of calls to a handful of candidates, combined with segregation and frequency filters ^{5 6}. You have no family segregation data on a single genome, so the gene list carries more of the load.

One warning about loss-of-function calls: a large fraction of predicted LoF variants in a healthy genome are annotation artifacts, in the last exon, on a non-expressed transcript, or in a gene with no LoF constraint. Check the gnomAD pLI / LOEUF for the gene before you spend time on the variant. LOEUF above 1.0 means the gene tolerates LoF and your “damaging” variant is probably inert.

7. Search each candidate by hand

Now the shortlist goes into a browser, one variant at a time. Four resources, each for something different.

ClinVar is the primary record. Search by the HGVS string from VEP (NM_000059.4:c.1234del) or by rs ID. Read the individual submissions. The aggregate label alone is not enough. Check submission dates and whether the submitters are clinical labs with assertion criteria.

ClinGen tells you whether the gene-disease relationship is validated at all, and whether a Variant Curation Expert Panel has a specification for this gene. The gene-level dosage sensitivity and validity scores are what keep you from chasing a “pathogenic” variant in a gene with disputed disease association.

VarSome aggregates roughly 30 databases behind one search box and runs an automated ACMG classifier on any variant you paste, including ones with no ClinVar record ⁷. URLs are predictable, so you can script them: https://varsome.com/variant/hg38/1-55039974-G-T. The value is breadth in one place. The caution is that the automated classification is a starting point. It is not a verdict. Published evaluation of the classifier in a diagnostic setting found that its rules need modification against gene-specific criteria before the output can be used clinically ⁸. Reviews of automated ACMG tools generally report substantial disagreement between tools and between tools and expert panels, particularly for variants of uncertain significance, which is precisely the category your shortlist is full of ⁹.

Franklin (Genoox) is a useful second automated classifier for the same reason a second compiler is useful: where two ACMG engines disagree, the variant needs human curation. Where both say benign and gnomAD says 2% and ClinVar has three two-star benign submissions, you can close the tab.

Record the following for each candidate.

• HGVS on the MANE transcript
• gnomAD grpmax AF and homozygote count
• ClinVar accession and review status
• gene-disease validity
• inheritance mode
• your zygosity

A heterozygous variant in a recessive gene means one thing, in a dominant gene something entirely different, and this is the distinction people get wrong most often when reading their own data.

8. Know where the line is

Nothing above is a diagnosis, and short-read WGS is not a clinical test. Coverage is uneven and indels in homopolymers are unreliable. Pseudogene-rich loci (PMS2, SMN1/2, CYP2D6 and HBA) are routinely miscalled by standard pipelines ^{2 3}. Sequencing data changes interpretation over time as databases are recurated, which means a shortlist is a snapshot rather than a final answer ¹⁰.

If any candidate survives step 7 and concerns you, the next step is a genetic counselor or clinical geneticist, and orthogonal confirmation in a CLIA/CAP-accredited lab. Do not act on a research-grade call, and do not make a health decision from a VarSome badge. Re-run steps 3 and 4 every six months against fresh ClinVar and gnomAD releases. The annotations move even when your genome does not.

Common problems

Zero ClinVar matches on a genome that clearly has some. Chromosome naming mismatch (chr1 vs 1), or you skipped normalization. Test by checking a known common variant with bcftools view -r chr1:55039974 norm.clinvar.vcf.gz.

Thousands of “pathogenic” hits. You matched on position only, or you are including CLNSIG=Conflicting_interpretations and zero-star submissions. Add the CLNREVSTAT filter.

Every candidate is in a gene you have never heard of, all olfactory receptors and HLA. Highly polymorphic and paralog-rich regions produce mapping artifacts. Exclude the ENCODE blacklist and segmental duplication regions with bcftools view -T ^blacklist.bed.

VEP consequences do not match the transcript your clinician uses. You picked a non-MANE transcript. Re-run with --mane --pick_allele_gene and report HGVS on the MANE Select accession.

Liftover output has half the expected variants. Wrong chain file direction, or the target FASTA does not match. Read the .unmap file that CrossMap writes.

Your CADD scores are all missing. The plugin file path is wrong or the tabix index is stale. VEP fails soft on plugins and will not error out.

Footnotes

1. Noam Hadar, Vadim Dolgin, Katya Oustinov, et al. VARista: a free web platform for streamlined whole-genome variant analysis across T2T, hg38, and hg19. Human Genetics, 2024. https://doi.org/10.1007/s00439-024-02671-4 ↩ ↩²

2. Neil Weisenfeld, Shuangye Yin, Ted Sharpe, et al. Comprehensive variation discovery in single human genomes. Nature Genetics, 2014. https://doi.org/10.1038/ng.3121 ↩ ↩²

3. Adrian V. Dalca, Michael Brudno. Genome variation discovery with high-throughput sequencing data. Briefings in Bioinformatics, 2010. https://doi.org/10.1093/bib/bbp058 ↩ ↩²

4. Alireza Tafazoli, John Mikros, Faeze Khaghani, et al. Pharmacovariome scanning using whole pharmacogene resequencing coupled with deep computational analysis and machine learning for clinical pharmacogenomics. Human Genomics, 2023. https://doi.org/10.1186/s40246-023-00508-1 ↩

5. Aayushi Srivastava, Diamanto Skopelitou, Beiping Miao, et al. Prioritization of predisposition genes for familial non-medullary thyroid cancer by whole-genome sequencing. European Journal of Endocrinology, 2025. https://doi.org/10.1093/ejendo/lvaf045 ↩

6. Aayushi Srivastava, Sara Giangiobbe, Abhishek Kumar, et al. Identification of Familial Hodgkin Lymphoma Predisposing Genes Using Whole Genome Sequencing. Frontiers in Bioengineering and Biotechnology, 2020. https://doi.org/10.3389/fbioe.2020.00179 ↩

7. Christos Kopanos, Vasilis Tsiolkas, Alexandros Kouris, et al. VarSome: the human genomic variant search engine. Bioinformatics, 2018. https://doi.org/10.1093/bioinformatics/bty897 ↩

8. Francesca Cristofoli, Elisa Sorrentino, Giulia Guerri, et al. Variant Selection and Interpretation: An Example of Modified VarSome Classifier of ACMG Guidelines in the Diagnostic Setting. Genes, 2021. https://doi.org/10.3390/genes12121885 ↩

9. Mireia Costa, Alberto García S., Ana León Palacio, et al. The promises and pitfalls of automated variant interpretation: a comprehensive review. Briefings in Bioinformatics, 2025. https://doi.org/10.1093/bib/bbaf545 ↩

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