How to Filter a Whole-Genome VCF Down to Variants Worth Reading

Oak
vcf variant-filtering wgs bcftools annotation
An iridescent winged reptile on a black pedestal, its glassy wing scales mostly dark with a few glowing amber, surrounded by hanging blank tags and lenses.

You will end up with three files: a high-confidence variant set (sample.pass.vcf.gz, roughly 3.5-4.5M variants for one human genome at 30x), an annotated table of those variants joined to gnomAD frequencies and predicted consequence, and a short candidate list, typically 20-80 rows, that you can read line by line. You need a single-sample VCF from a 30x WGS run (ideally a gVCF you genotyped yourself, or the joint-called VCF your provider returned), the matching reference FASTA, about 60 GB of disk, and bcftools ≥1.17, slivar ≥0.3, and either ANNOVAR or VEP. A laptop with 16 GB of RAM is enough. Everything below is measurement and interpretation. It does not tell you what any variant means for your health, and any variant you plan to act on needs confirmation in a clinical lab and a conversation with a genetic counselor or physician.

1. Look at what you have before you filter anything

Do not run a single filter until you know the shape of the file. Header first:

bcftools view -h sample.vcf.gz | grep -E '^##(FILTER|INFO|FORMAT|source|reference)' | head -60

You are checking four things. Was this called by GATK HaplotypeCaller, DeepVariant, or DRAGEN? Are the FILTER values real (PASS, VQSRTrancheSNP99.00to99.90, RefCall) or is every record .? Does the INFO field carry VQSLOD (GATK VQSR), or GQ/QUAL only (DeepVariant, which does not emit VQSR annotations and needs different thresholds)? And is the reference build hg19 or hg38? Build mismatch between your VCF and your annotation database is the most common way to produce a silently wrong answer, because coordinates shift by a few hundred to a few thousand bases in most regions and you will get annotations for the wrong gene without any error message.

Then the basic counts:

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

Sanity ranges for one 30x human genome, GRCh38, single sample: 4.0-5.0M SNVs, 500-900k indels, Ti/Tv ratio 2.0-2.1 genome-wide (2.8-3.0 if you subset to exons). A Ti/Tv below 1.8 means you have a pile of false-positive transversions and your downstream candidate list will be mostly noise. A total variant count above 6M usually means a contaminated or mixed sample, or that you are looking at a multi-sample file.

bcftools query -l sample.vcf.gz | wc -l   # confirm it is one sample

2. Apply site-level and genotype-level hard filters

Two filtering layers exist and people conflate them. Site-level filters ask whether a variant exists at this position at all, using INFO annotations. Genotype-level filters ask whether this sample’s specific call is trustworthy, using FORMAT fields. You need both. The canonical filter set in the literature is built on genotype quality, depth, allele balance, and population frequency 1.

If your VCF came from GATK with VQSR already applied, respect it and keep only the tranches you want:

bcftools view -f PASS -Oz -o sample.site.vcf.gz sample.vcf.gz
bcftools index -t sample.site.vcf.gz

If FILTER is empty (raw HaplotypeCaller output, or a provider who skipped VQSR), hard-filter. SNVs and indels need different thresholds because indel calling is noisier:

# SNVs
bcftools view -v snps sample.vcf.gz \
| bcftools filter -s LowQual -m + \
    -e 'QUAL<30 || INFO/QD<2.0 || INFO/FS>60.0 || INFO/MQ<40.0 || INFO/MQRankSum<-12.5 || INFO/ReadPosRankSum<-8.0' \
    -Oz -o snv.filt.vcf.gz

# indels
bcftools view -v indels,mnps,other sample.vcf.gz \
| bcftools filter -s LowQual -m + \
    -e 'QUAL<30 || INFO/QD<2.0 || INFO/FS>200.0 || INFO/ReadPosRankSum<-20.0' \
    -Oz -o indel.filt.vcf.gz

bcftools index -t snv.filt.vcf.gz; bcftools index -t indel.filt.vcf.gz
bcftools concat -a -Oz -o sample.site.vcf.gz snv.filt.vcf.gz indel.filt.vcf.gz
bcftools index -t sample.site.vcf.gz

Those numbers are GATK’s published hard-filter recommendations, and they are deliberately loose: they are tuned to throw away almost nothing real. Do not tighten them by guesswork. If you want a tighter set, use VQSR or a machine-learned filter rather than moving QD to 5.

Now the genotype layer. This is where most of the junk in a candidate list comes from:

bcftools view -f PASS sample.site.vcf.gz \
| bcftools filter -i 'FMT/GQ>=20 && FMT/DP>=10 && FMT/DP<=250' \
| bcftools view -e 'GT="0/0" || GT="./."' \
    -Oz -o sample.pass.vcf.gz
bcftools index -t sample.pass.vcf.gz

GQ ≥ 20 means the called genotype is about 100x more likely than the next best. DP ≥ 10 is the floor below which heterozygous calls start dropping out from sampling alone: at 10x, the chance of seeing zero alt reads at a true het is about 0.1%, but the chance of seeing 2 or fewer is a few percent, which is why low-depth hets are unreliable. The DP ≤ 250 cap removes collapsed repeats and segmental duplications, where multiple genomic copies pile onto one locus and produce confident-looking, entirely false heterozygous calls.

Allele balance is the single most useful genotype filter and most people skip it. A true het should have alt fraction near 0.5:

bcftools filter -i '(FMT/AD[0:1]) / (FMT/AD[0:0] + FMT/AD[0:1]) > 0.22' \
  sample.pass.vcf.gz -Oz -o sample.ab.vcf.gz

For homozygous-alt calls, require alt fraction above 0.9. Calls between 0.05 and 0.20 alt fraction are mostly mapping artifacts, index hopping, or somatic mosaicism, and you cannot tell which from the VCF alone.

Expect to lose 3-8% of variants across steps 2. If you are losing 30%, something upstream is wrong: check for a truncated BAM, an unindexed reference, or a sample that was sequenced at 15x and labeled 30x.

3. Restrict to regions where short reads work

About 2-3% of the genome is inaccessible to 150 bp paired-end reads, and that 2-3% generates a wildly disproportionate share of false positives. Intersect with a mappability or accessibility mask. The GA4GH/GIAB high-confidence regions BED is the one we use:

bcftools view -R GRCh38_notinalldifficultregions.bed.gz \
  sample.ab.vcf.gz -Oz -o sample.acc.vcf.gz
bcftools index -t sample.acc.vcf.gz

This removes homopolymers longer than 10 bp, tandem repeats, segmental duplications, and the MHC. You lose real variation in those regions, which matters if you care about HLA typing or specific repeat expansions, and those need dedicated callers anyway (ExpansionHunter for repeats, HLA-LA or OptiType for HLA). For general prioritization, cutting them is the right trade.

4. Annotate with consequence and population frequency

Annotation is where filtering becomes interpretation. You need, at minimum: gene and transcript, predicted consequence, gnomAD allele frequency, and ClinVar classification. ANNOVAR’s table-based approach is fast and well-documented 2, and its table_annovar.pl with multiple protocols in one pass is the least painful route:

bcftools norm -m-both -f GRCh38.fa sample.acc.vcf.gz -Oz -o sample.norm.vcf.gz

table_annovar.pl sample.norm.vcf.gz humandb/ \
  -buildver hg38 \
  -out sample \
  -remove \
  -protocol refGene,gnomad41_genome,clinvar_20240917,dbnsfp47a \
  -operation g,f,f,f \
  -nastring . \
  -vcfinput \
  -polish \
  -thread 8

bcftools norm -m-both is not optional. Multi-allelic records annotate incorrectly in almost every tool because the annotation gets attached to the first ALT and copied. -polish fixes ANNOVAR’s indel left-alignment against the reference, which changes HGVS notation for a meaningful fraction of indels.

VEP is the alternative and gives better transcript-level detail plus a plugin ecosystem (SpliceAI, CADD, AlphaMissense). It is slower and needs a 25 GB cache. If you only want a web-based pass over a few hundred variants, SNPnexus handles annotation against multiple gene and regulatory datasets without a local install 3. For a full genome, run locally.

Output is sample.hg38_multianno.txt, a tab-delimited table with one row per variant and 100+ columns. This is the file you will work in.

5. Reduce to candidates with frequency and consequence filters

Now the ordering that matters. Filter on frequency first, because it is the highest-yield single step, then on consequence.

For a rare-disease framing, the standard cut is gnomAD popmax allele frequency below 0.1% for dominant conditions and below 1% for recessive 4. For a healthy-individual baseline, we use 0.5% popmax as a first pass and look at what survives.

Using slivar, which is built exactly for this and is fast:

slivar expr \
  --vcf sample.norm.vcf.gz \
  --gnotate gnomad.v4.1.hg38.zip \
  --info 'INFO.gnomad_popmax_af < 0.005 && variant.FILTER == "PASS"' \
  --sample-expr 'het:sample.het && sample.GQ > 20 && sample.DP > 10 && sample.AB > 0.25 && sample.AB < 0.75' \
  --sample-expr 'hom:sample.hom_alt && sample.GQ > 20 && sample.DP > 10 && sample.AB > 0.90' \
  --pass-only \
  -o sample.rare.vcf.gz

Then consequence. Keep only variants predicted to change protein or splicing:

awk -F'\t' 'NR==1 || $9 ~ /^(frameshift|stopgain|stoploss|startloss|nonsynonymous|splicing)/' \
  sample.hg38_multianno.txt > sample.coding.txt

(Column 9 is ExonicFunc.refGene in the layout above. Verify with head -1 | tr '\t' '\n' | nl. Column positions shift between ANNOVAR versions and protocol lists.)

A 30x genome typically leaves 300-800 rare coding variants after this. That is still too many to read, so add two more constraints. First, restrict to genes with established disease association rather than the whole genome: a ClinGen-curated gene list, or the ~4,500 genes in a clinical exome panel. Gene-level restriction is the filter that does the most to reduce false leads 5. Second, prefer variants with supporting evidence: a ClinVar entry with two-star review status or better, a CADD phred above 20, or a SpliceAI delta score above 0.5 for splice candidates.

# join annotated table to a gene list
awk -F'\t' 'NR==FNR{g[$1]; next} FNR==1 || ($7 in g)' \
  clingen_definitive_genes.txt sample.coding.txt > sample.candidates.txt
wc -l sample.candidates.txt

That should land you between 20 and 80 rows. Read every one.

6. Verify the survivors in the pileup

Every variant you take seriously gets looked at. Load the BAM in IGV, jump to the locus, and check: do alt reads appear on both strands, at varied positions within the read, with good mapping quality? Artifacts cluster at read ends, appear on one strand only, or sit inside a soft-clipped tail. Three to five minutes per variant, and it will kill a third of a typical candidate list.

samtools tview -p chr17:43094000 -d T sample.bam GRCh38.fa

or better, samtools mpileup -r chr17:43094000-43094100 -f GRCh38.fa sample.bam and read the column.

Research-grade WGS is not a clinical result. If a variant looks real and consequential, the next step is orthogonal confirmation (Sanger or a clinical panel) in a CAP/CLIA lab, ordered through a clinician. Do not act on a research VCF.

Common problems

Most of your survivors will be VUS, and you cannot resolve them yourself. In one large hereditary-disease testing series, variants of uncertain significance were reported at rates that dwarf pathogenic findings, and the rate varies substantially by gene and by ancestry 6. A VUS means the lab has insufficient evidence to place the variant on either side, not that the variant is halfway pathogenic. The correct default is to treat it as uninformative until evidence accumulates 7. Reclassification happens: many VUS are eventually moved to benign, a small minority to pathogenic, and the direction is not predictable from the variant’s appearance 8. Large-scale functional assays (multiplexed assays of variant effect) plus better population data are expected to shrink the VUS pool considerably over the next several years, but the problem will not disappear 9.

A variant can be benign for straightforward mechanical reasons. It sits in a gene where loss of function has no phenotype. It changes an amino acid to one with the same chemistry at a position that is not conserved. It is a synonymous change far from a splice boundary. Or it is common: a variant at 5% frequency in any population is almost certainly not causing a severe dominant condition, because the condition would then be far more prevalent than it is. Frequency-based exclusion is a strong argument, which is why it goes first in the pipeline.

Ancestry bias in the frequency databases. gnomAD is still dominated by European-ancestry samples. If your ancestry is underrepresented, benign variants common in your population may show low global frequency and survive your filter. This is a documented driver of higher VUS rates in non-European patients 6. Use popmax rather than global AF, and treat “not in gnomAD” as weaker evidence than it looks.

Over-filtering hides the answer. Every threshold you tighten trades sensitivity for specificity, and the causal variant is sometimes the one with DP 9 or GQ 18. A useful practice: keep the intermediate files, and when your candidate list is empty, relax one filter at a time and note what appears. Published filtering guides emphasize this iterative loosening rather than a single fixed cascade 4.

Multi-allelic and non-normalized records. If you skipped bcftools norm -m-both -f ref.fa, your frequency lookups will silently miss. A left-aligned indel in gnomAD will not match a right-aligned indel in your VCF, so it will look novel and rare. Always normalize before annotating and before any database join.

Structural variants and repeat expansions are absent. Short-read SNV/indel callers do not report deletions above ~50 bp, inversions, or expanded repeats. If nothing in your candidate list explains a phenotype you care about, that is a reason to run manta or smoove for SVs and ExpansionHunter for known repeat loci, not a reason to conclude nothing is there.

Confusing prediction scores with evidence. CADD, REVEL, and AlphaMissense outputs are priors, not findings. A CADD of 32 on a variant in a gene with no disease association is not a result. In ACMG-style frameworks, computational evidence carries supporting weight at most, and stacking three correlated predictors does not add three pieces of evidence. The reporting and communication burden around uncertain findings is real and well-described, and it is a reason to bring a clinician in before you go far down any interpretive path 10.

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

  1. Joke Reumers, Peter De Rijk, Hui Zhao, et al. Optimized filtering reduces the error rate in detecting genomic variants by short-read sequencing. Nature Biotechnology, 2011. https://doi.org/10.1038/nbt.2053 ↩

  2. Hui Yang, Kai Wang. Genomic variant annotation and prioritization with ANNOVAR and wANNOVAR. Nature Protocols, 2015. https://doi.org/10.1038/nprot.2015.105 ↩

  3. Jorge Oscanoa, Lavanya Sivapalan, Emanuela Gadaleta, et al. SNPnexus: a web server for functional annotation of human genome sequence variation (2020 update). Nucleic Acids Research, 2020. https://doi.org/10.1093/nar/gkaa420 ↩

  4. Mahjoubeh Jalali Sefid Dashti, Junaid Gamieldien. A Practical Guide To Filtering and Prioritizing Genetic Variants. BioTechniques, 2017. https://doi.org/10.2144/000114492 ↩ ↩2

  5. David Salgado, Matthew I. Bellgard, Jean-Pierre Desvignes, et al. How to Identify Pathogenic Mutations among All Those Variations: Variant Annotation and Filtration in the Genome Sequencing Era. Human Mutation, 2016. https://doi.org/10.1002/humu.23110 ↩

  6. Elaine Chen, Flavia M. Facio, Kerry W. Aradhya, et al. Rates and Classification of Variants of Uncertain Significance in Hereditary Disease Genetic Testing. JAMA Network Open, 2023. https://doi.org/10.1001/jamanetworkopen.2023.39571 ↩ ↩2

  7. Karen E. Weck. Interpretation of genomic sequencing: variants should be considered uncertain until proven guilty. Genetics in Medicine, 2018. https://doi.org/10.1038/gim.2017.269 ↩

  8. Giulia Federici, Silvia Soddu. Variants of uncertain significance in the era of high-throughput genome sequencing: a lesson from breast and ovary cancers. Journal of Experimental & Clinical Cancer Research, 2020. https://doi.org/10.1186/s13046-020-01554-6 ↩

  9. Douglas M. Fowler, Heidi L. Rehm. Will variants of uncertain significance still exist in 2030?. The American Journal of Human Genetics, 2024. https://doi.org/10.1016/j.ajhg.2023.11.005 ↩

  10. Lily Hoffman-Andrews. The known unknown: the challenges of genetic variants of uncertain significance in clinical practice. Journal of Law and the Biosciences, 2017. https://doi.org/10.1093/jlb/lsx038 ↩