How to Use Allele Frequency Databases With Your Own Genome

By the end of this guide you will have a normalized, frequency-annotated VCF of your own genome in which every variant carries population allele frequencies from a reference cohort. You will also have a filtered subset of variants that are rare or absent in every reference population. Finally you will have a per-variant Hardy-Weinberg statistic you can use to flag calls that are probably artifacts rather than biology. You need a germline VCF from whole-genome or whole-exome sequencing, ideally with genotype quality and depth fields intact. You also need the reference FASTA the calls were made against and about 100 GB of free disk. The tools are bcftools 1.19 or later, echtvar, and plink2. Everything below runs on a laptop except the initial download of the gnomAD frequency tables, which is bandwidth-bound.

1. Normalize your VCF before you touch a database

Frequency annotation is a join on a key made of chromosome, position, reference allele and alternate allele. If your representation of an indel differs from the database’s by a single base of left-padding, the join silently fails and the variant looks novel. This is the single most common reason a personal-genome pipeline reports thousands of “never-before-seen” variants that are in fact common.

Split multiallelic records, left-align indels against the reference, and drop records whose REF does not match the reference sequence:

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

-m -any splits A > C,G into two records so each has one ALT allele. --check-ref w writes a warning and fixes the REF where it can, rather than exiting. -d exact removes duplicate records introduced by the split. Then confirm your build: the database you annotate with must be on the same assembly. A quick check is the length of chr1 in the header, 248,956,422 for GRCh38 and 249,250,621 for GRCh37. Mixing builds produces frequencies that are wrong rather than missing, which is worse.

Also confirm that your VCF distinguishes homozygous-reference sites from uncalled sites. A single-sample VCF that contains only variant lines cannot tell you whether a position was covered at 40x and matched the reference or was never sequenced at all. If you have a gVCF, keep it: reference confidence blocks are what let you say “I do not carry this variant” as opposed to “I have no data here.” The distinction matters when a downstream question depends on absence.

2. Choose the reference cohort deliberately

The frequency you get depends entirely on who was sequenced. The main general-purpose options, and when we would reach for each:

gnomAD is the default. Version 4 aggregates roughly 800,000 exomes and genomes, reports frequencies split by genetic ancestry group, and publishes allele number (AN) per site so you can see how well a position was covered across the cohort. Use the genome subset for non-coding regions and the joint frequencies for coding variants. Its weakness is that it is an aggregate of case-control and biobank studies, so it is not a random sample of humanity. Severe pediatric disease alleles are depleted in it but not absent.

The 1000 Genomes Project phase 3 remains useful because it gives you individual-level genotypes in 26 populations rather than summary counts ¹. That means you can compute your own frequencies, subsets, and haplotypes. Its 2,504 samples are too few for anything rarer than about 0.1%, and the original low-coverage design means genotypes at rare sites were imputed rather than directly observed. The 30x re-sequencing of the same panel largely fixes that.

Population-specific databases matter when your ancestry is underrepresented in gnomAD. The Human Genetic Variation Database reports variant frequencies derived from Japanese exomes and is the reference we would check for a variant that looks rare globally but might be common in East Asia ². The Thousand Polish Genomes database plays the same role for Central and Eastern European ancestry, and its authors found a substantial number of variants at appreciable frequency in Poland that were rare or absent in pan-European aggregates ³. ALFRED collected allele frequencies across a large number of anthropologically defined populations ⁴. It covers older and more finely partitioned population samples, which makes it useful for teaching and for classical markers.

For HLA and other immune genes, general-purpose VCF databases are the wrong tool, and step 7 covers the alternative.

3. Annotate with echtvar rather than a VCF-to-VCF join

You can annotate with bcftools annotate -a gnomad.vcf.bgz, and for a single region that is fine. For a whole genome it is slow, because the annotation source is a multi-gigabyte compressed VCF that must be parsed. echtvar pre-encodes the frequency fields into a compact archive keyed by position and allele, and annotates a whole genome in a couple of minutes.

Build the archive once. Write a small JSON config describing the fields you want to carry over:

[
  {"field": "AF", "alias": "gnomad_AF", "multiplier": 2000000, "missing_value": -1},
  {"field": "AF_grpmax", "alias": "gnomad_AF_grpmax", "multiplier": 2000000, "missing_value": -1},
  {"field": "AN", "alias": "gnomad_AN"},
  {"field": "nhomalt", "alias": "gnomad_nhomalt"}
]

The multiplier is how echtvar stores a float as an integer. At 2,000,000 the smallest representable nonzero frequency is 5e-7, finer than the resolution of an 800,000-person cohort, so nothing meaningful is lost. Then:

echtvar encode gnomad_v4_af.zip conf.json gnomad.joint.v4.1.sites.chr*.vcf.bgz
echtvar anno -e gnomad_v4_af.zip mine.norm.vcf.gz mine.gnomad.bcf

Check that the join worked before you trust it. The fraction of your variants with a nonzero gnomAD frequency should be high. For a single human genome against gnomAD v4, we expect the large majority of single-nucleotide variants to have a frequency record. The unmatched remainder is concentrated in indels, repeats, and low-complexity regions. If more than a few percent of your SNVs are unmatched, you have a build or normalization mismatch. A genome full of private variation is not the explanation.

Now filter. To pull variants that are rare everywhere, use the maximum frequency across ancestry groups rather than the global average, because a variant at 4% in one group and absent elsewhere has a low global frequency and is not rare:

bcftools view -i 'gnomad_AF_grpmax < 0.001 || gnomad_AF_grpmax = "."' \
  -Ob -o mine.rare.bcf mine.gnomad.bcf

Also require that the site was well covered in the reference cohort. A variant with gnomad_AN far below the cohort maximum sits in a region where the database itself has poor data, and its apparent rarity is an artifact of missing calls rather than a statement about the population.

4. Compute allele frequencies yourself

If you have genotypes for a cohort, the frequency is a count, and the formula is worth stating plainly. Consider a biallelic site with genotype counts N_AA, N_Aa, and N_aa in N diploid individuals. The frequency of allele A is

p = (2·N_AA + N_Aa) / (2·N)

and q = 1 − p. The denominator is the allele number: two per individual at an autosomal site, so 2N in total, which is what VCF files report as AN. The numerator is the allele count, AC. So AF = AC/AN, and the “most common allele” at a site is simply the one with the largest count, which is not always the reference allele. Sites where the alternate allele is the major one are common enough that any pipeline filtering on “is it reference” will misbehave.

With bcftools you can fill these in for any multi-sample VCF:

bcftools +fill-tags cohort.vcf.gz -Oz -o cohort.tags.vcf.gz -- -t AN,AC,AF,NS,HWE,ExcHet

A frequency of 1.0 means every allele observed at that site was the alternate allele: the variant is fixed in the sample, and no reference allele was seen. In a personal genome this usually means the reference assembly carries the minor allele at that position, and you and everyone else in the cohort carry the common one. It is not a sign of anything unusual about you. Frequencies of exactly 0.0 in a database can likewise mean either “tested and never seen” or “no data,” which is why AN matters more than AF alone.

5. Test Hardy-Weinberg equilibrium and know what it tells you

Hardy-Weinberg equilibrium describes the genotype frequencies expected from allele frequencies alone under random mating: p² homozygous for A, 2pq heterozygous, q² homozygous for a. You test it by comparing observed genotype counts to those expectations. Use a chi-square test on one degree of freedom, or preferably an exact test. The exact test is what you want for rare variants where the expected homozygote count is below five. The classical assumptions are five. There is no selection at the locus, no mutation, and no migration. The population is effectively infinite so that drift is negligible, and mating is random with respect to the genotype. Two further conditions are usually left implicit, namely that the organism is diploid with non-overlapping generations and that the locus is autosomal, which is why you exclude the X chromosome outside the pseudoautosomal regions from any genome-wide HWE scan.

In practice you are almost never testing HWE to learn population genetics. You are using it as a quality-control filter, because sequencing and genotyping artifacts violate it in characteristic ways. Mapping errors between paralogous regions inflate heterozygosity. Allelic dropout and reference bias inflate homozygosity. Genome-wide studies with next-generation sequence data show that HWE deviations cluster in exactly the regions you would suspect, and that the test flags systematic genotyping error rather than interesting biology at most sites ⁵. Large-scale sequencing databases show the same pattern, which is part of why HWE filtering has become routine in variant-database curation ⁶.

To run it:

plink2 --vcf cohort.vcf.gz --hwe 1e-10 midp keep-fewhet \
  --autosome --geno 0.05 --maf 0.001 \
  --make-bed --out cohort.qc
plink2 --bfile cohort.qc --hardy midp --out cohort.hwe

The midp option applies the mid-p adjustment to the exact test, which is less conservative at low counts. keep-fewhet retains sites that deviate because they have too few heterozygotes, on the argument that heterozygote deficiency in a case-control cohort can reflect real biology while heterozygote excess is almost always a mapping artifact. A threshold of 1e-10 rather than the often-quoted 1e-6 is deliberate: with millions of variants tested, a lenient threshold discards real sites, and you are looking for gross failures.

Two cautions. First, population structure breaks the test: a cohort that mixes ancestries shows heterozygote deficiency at any site whose frequency differs between the groups, an effect that has nothing to do with genotyping quality. Your cohort may be structured, or your data may be low-depth so that genotypes are uncertain. In either case use a method that models both. Such a method tests HWE from genotype likelihoods with admixture proportions estimated from the data rather than from called genotypes ⁷. Second, HWE is a property of a population sample. It does not apply to an individual. You cannot test your own genome for HWE. What you can do is inherit the test: filter your variants to those that pass HWE in the reference database, which is a cheap and effective way to remove recurrent artifact sites.

6. Interpret frequency without overreaching

Once annotated, the temptation is to read rarity as consequence. Resist it. Every healthy person carries on the order of tens of thousands of variants absent from any database and a few hundred predicted loss-of-function variants, most of them with no discernible effect. Frequency is a prior, not a conclusion: a variant at 5% in your ancestry group is very unlikely to cause a severe early-onset condition, while a variant seen once in 800,000 people is merely unexplained.

The useful move is to combine frequency with independent evidence. Ask whether the gene is known to be associated with a phenotype and whether the variant class is one that plausibly disrupts the protein. Then ask whether it segregates in your family and whether your own RNA and protein measurements show anything at that locus. And when a specific variant looks clinically meaningful, the next step is a clinical-grade confirmatory assay ordered and interpreted by a clinician or certified genetics laboratory. Research-grade sequencing and public frequency databases are not a diagnostic and should not be treated as one. Note also that reported frequencies shift between database versions as cohorts grow and calling pipelines change, so record the version you used alongside any result you plan to revisit ⁸.

7. HLA and immune loci need a different database

The HLA region is the most polymorphic part of the human genome. Its alleles are defined as full-length sequence haplotypes with four-field names such as A*02:01:01:01. They are not defined as individual positions in a VCF. Standard variant calling against the linear reference performs poorly there, so you type HLA separately. Tools for this include HLA-LA, OptiType, and arcasHLA. The last of these types from RNA-seq reads and is convenient if you already have a transcriptome.

Once you have allele calls, the reference for their frequency is the Allele Frequency Net Database ⁹. It curates allele, haplotype, and genotype frequencies across several thousand population samples. The coverage includes HLA and other immune gene systems such as KIR and cytokine loci. Query your specific alleles against populations that match your ancestry, because HLA allele frequencies vary more sharply across populations than almost any other locus. Treat any single typing call with appropriate suspicion until confirmed by a second method, since ambiguity within an allele group is common.

Common problems

A large fraction of variants come back with no frequency annotation. Almost always a build or normalization problem. Confirm the assembly on both sides, rerun bcftools norm with -f and --check-ref w, and check whether your contig naming uses chr1 or 1 consistently with the annotation source.

Common variants show up as rare. Check that you are filtering on the ancestry-group maximum rather than the global frequency, and that you have not restricted the annotation to the exome subset when your variant is intronic. A coding-only frequency source gives AN near zero outside exons, and near-zero AN looks like near-zero frequency.

Your rare-variant count is implausibly high. Look at the distribution across the genome. If the excess concentrates in segmental duplications, centromeric regions, and low-complexity sequence, it is mapping error. Intersect with a high-confidence region mask and re-count.

HWE fails everywhere in a small cohort. With fewer than a few hundred samples the exact test has little power at low frequency, and the p-value distribution is dominated by discreteness. Restrict the scan to sites with minor allele frequency above about 1% before drawing conclusions, and remember that a structured sample produces heterozygote deficiency for reasons unrelated to data quality ⁷.

A variant is absent from gnomAD but present in a population-specific database. This is expected and informative rather than contradictory. Underrepresented ancestries carry variants that pan-cohort aggregates miss, which is the whole reason the population-specific resources exist ^{2 3}. When the two disagree, weight the database whose sample matches your ancestry.

Indels disagree between your calls and the database. Left-alignment handles most cases, but complex indels and multi-nucleotide substitutions can be represented in genuinely different ways that no normalization reconciles. For a small number of variants you care about, inspect the alignments directly in IGV and compare against the database’s own representation rather than trusting the join.

Footnotes

1. The 1000 Genomes Project Consortium. A map of human genome variation from population-scale sequencing. Nature, 2010. https://doi.org/10.1038/nature09534 ↩

2. Koichiro Higasa, Noriko Miyake, Jun Yoshimura, et al. Human genetic variation database, a reference database of genetic variations in the Japanese population. Journal of Human Genetics, 2016. https://doi.org/10.1038/jhg.2016.12 ↩ ↩²

3. Elżbieta Kaja, Adrian Lejman, Dawid Sielski, et al. The Thousand Polish Genomes—A Database of Polish Variant Allele Frequencies. International Journal of Molecular Sciences, 2022. https://doi.org/10.3390/ijms23094532 ↩ ↩²

4. H. Rajeevan, U. Soundararajan, J. R. Kidd, et al. ALFRED: an allele frequency resource for research and teaching. Nucleic Acids Research, 2011. https://doi.org/10.1093/nar/gkr924 ↩

5. Jan Graffelman, Deepti Jain, Bruce Weir. A genome-wide study of Hardy–Weinberg equilibrium with next generation sequence data. Human Genetics, 2017. https://doi.org/10.1007/s00439-017-1786-7 ↩

6. Nikita Abramovs, Andrew Brass, May Tassabehji. Hardy-Weinberg Equilibrium in the Large Scale Genomic Sequencing Era. Frontiers in Genetics, 2020. https://doi.org/10.3389/fgene.2020.00210 ↩

7. Jonas Meisner, Anders Albrechtsen. Testing for Hardy–Weinberg equilibrium in structured populations using genotype or low‐depth next generation sequencing data. Molecular Ecology Resources, 2019. https://doi.org/10.1111/1755-0998.13019 ↩ ↩²

8. André Altmann, Peter Weber, Daniel Bader, et al. A beginners guide to SNP calling from high-throughput DNA-sequencing data. Human Genetics, 2012. https://doi.org/10.1007/s00439-012-1213-z ↩

9. Faviel F. Gonzalez-Galarza, Antony McCabe, Eduardo J. Melo dos Santos, et al. Allele Frequency Net Database. Methods in Molecular Biology, 2024. https://doi.org/10.1007/978-1-0716-3874-3_2 ↩