Compare my VCF to gnomAD variants

Question

I have a VCF with small variant calls against HG38 and I would like to determine which of those calls are present in the gnomAD database.

Is there an existing tool that can do this? Should I be looking at variant annotation tools? or is this something that I would look at BCFtools for?

Answer 1

You can use vcfanno to annotate your vcf with the gnomad vcf file. There's an example config file for annotating with gnomad (along with other databases such as dbSNP) here.

Answer 2

I haven't checked how fast this is but you could try using Hail (disclaimer: I'm a maintainer). You'll need to install the GCS connector first to read from Google Cloud Storage:

curl https://broad.io/install-gcs-connector | python3

import hail as hl

mt = hl.import_vcf('your.vcf')
mt = hl.split_multi_hts(mt)  # use the same biallelic representation that gnomAD uses
mt.write('data.mt')  # store in Hail's native format for faster queries
mt = hl.read_matrix_table('data.mt')
sites = hl.read_table(
    'gs://gcp-public-data--gnomad/release/3.1.2/ht/genomes/gnomad.genomes.v3.1.2.sites.ht'
)
mt = mt.anti_join_rows(sites)

That'll filter down to just those variants (defined by chromosome, position, and alleles) which are not present in gnomAD v3.

Answer 3

If you have sufficient available disk space, you could download the files from the gnomad downloads page and use a bedtools intersect command to find the overlapping variants.

However, the recommended platforms for gnomAD are to use hail utilities developed by them and others, due to the constraints of working with ~1TB of data.

I am not an expert but potentially someone else knows more about using the site hail data.