Gnomad WGS AFs are often ~10x the Gnomad WES AF

[stu2]

Hi all,

First, apologies for the cross-posting - I’ve posted this in Biostars but couldn’t get to the bottom of it.

I am trying to investigate population allele frequencies from Gnomad v4.1 but I have found that I can’t get the population AFs from the 76,215 full-genome sequences to roughly agree with the AFs from the 730,947 WES samples. There seems to be a sub-set for which the reported AF from WGS is ~10x that from the WES (where both are reported). It’s only variants with WES population AF < 0.01 that seem to be affected (see scatter plot)

Here is my code, looking only at chromosome 22 (because it is small):

1. downloaded gnomad.genomes.v4.1.1.sites.chr22.20M_filt.vcf.gz and gnomad.exomes.v4.1.1.sites.chr22.vcf.bgz (with their .tbi files) from gnomad data then:

2. Tabulate AFs, only look at region in Chr22: 10MB → 30MB

bcftools view --apply-filters '.,PASS' -e 'lcr=1 || variant_type!="snv" ' -r chr22:10000000-30000000 -O z -o  gnomad.exomes.v4.1.1.sites.chr22.10M_30M_filt.vcf.gz   ../gnomad.exomes.v4.1.1.sites.chr22.vcf.bgz &
bcftools view --apply-filters '.,PASS' -e 'lcr=1 || variant_type!="snv" ' -r chr22:10000000-30000000 -O z -o  gnomad.genomes.v4.1.1.sites.chr22.10M_30M_filt.vcf.gz   ../gnomad.genomes.v4.1.1.sites.chr22.vcf.bgz

# rename AF to WES_AF in exome calls
echo '##INFO=<ID=WES_AF,Number=A,Type=Float,Description="gnomAD allele frequency from WES">' > AF_header_exome.txt
echo 'INFO/AF INFO/WES_AF' > rename_annots_exome.txt
bcftools annotate -h AF_header_exome.txt --rename-annots rename_annots_exome.txt -Oz -o gnomad.exomes.v4.1.1.sites.chr22.10M_30M_filt_WES_AF.vcf.gz  gnomad.exomes.v4.1.1.sites.chr22.10M_30M_filt.vcf.gz
tabix gnomad.exomes.v4.1.1.sites.chr22.10M_30M_filt_WES_AF.vcf.gz
tabix gnomad.genomes.v4.1.1.sites.chr22.10M_30M_filt.vcf.gz

# now transfer AF to exome:
bcftools annotate -a gnomad.genomes.v4.1.1.sites.chr22.10M_30M_filt.vcf.gz  -c INFO/AF -Oz -o gnomad.exomes.v4.1.1.sites.chr22.10M_30M_filt_WGSanno.vcf.gz  gnomad.exomes.v4.1.1.sites.chr22.10M_30M_filt_WES_AF.vcf.gz
tabix gnomad.exomes.v4.1.1.sites.chr22.10M_30M_filt_WGSanno.vcf.gz
bcftools query -H  -f '%CHROM\t%POS\t%REF\t%ALT\t%WES_AF\t%AF\t%AS_culprit\t%AS_VarDP\t%VarDP\t%MQRankSum\t%MQ\t%outside_ukb_capture_region\t%outside_broad_capture_region\n' gnomad.exomes.v4.1.1.sites.chr22.10M_30M_filt_WGSanno.vcf.gz > testtab4.tsv

3. Plot in R:

chr22 <- read.table('testtab4.tsv', na.strings = c('.'),
                    header = T, 
                    comment.char = '', 
                    sep='\t',
                    colClasses = c('character', 'numeric', 'character','character', 'numeric','numeric','character', 'numeric','numeric','numeric','numeric','character','character'))
colnames(chr22) <- c("CHROM" , "POS","REF", "ALT", "WES_AF", "WGS_AF", "AS_culprit", "AS_VarDP" , "VarDP", "MQRankSum", "MQ", "outside_ukb_capture_region", "outside_broad_capture_region")
chr22$WGS_WES_AF_RAT <- log10(chr22$WGS_AF/chr22$WES_AF)
MIN_AF <- 0.00001
# only look at variants detectable in both WGS and WES:
chr22_filt <- chr22[chr22$WES_AF > MIN_AF & chr22$WGS_AF > MIN_AF,]
hist(chr22_filt$WGS_WES_AF_RAT,200)
ggplot2::ggplot(chr22_filt) + 
  geom_point(aes(x=log10(WES_AF), y=log10(WGS_AF), col = outside_broad_capture_region), size=0.5) + 
  theme_bw()