Hi all,

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).

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:

# subset to reasonable size; filter out AC=0 entries
bcftools view gnomad.genomes.v4.1.1.sites.chr22.vcf.bgz -r chr22:10000000-20000000  | grep -v 'AC=0;AN=' |  bgzip --threads 8  >  gnomad.genomes.v4.1.1.sites.chr22.20M_filt.vcf.gz
bcftools view gnomad.exomes.v4.1.1.sites.chr22.vcf.bgz  -r chr22:10000000-20000000  | grep -v 'AC=0;AN=' |  bgzip --threads 8  >  gnomad.exomes.v4.1.1.sites.chr22.20M_filt.vcf.gz
# rename AF to gnomADg_AF in genome calls
echo '##INFO=<ID=gnomADg_AF,Number=A,Type=Float,Description="gnomAD allele frequency from WGS">' > AF_header.txt
echo 'INFO/AF INFO/gnomADg_AF' > rename_annots.txt
bcftools annotate -h AF_header.txt --rename-annots rename_annots.txt -Oz -o gnomad.genomes.v4.1.1.sites.chr22.20M_filt_gnomADg_AF_anno.vcf.gz  gnomad.genomes.v4.1.1.sites.chr22.20M_filt.vcf.gz
tabix gnomad.genomes.v4.1.1.sites.chr22.20M_filt_gnomADg_AF_anno.vcf.gz
# now transfer gnomADg_AF field to exome where variant matches:
bcftools annotate -a gnomad.genomes.v4.1.1.sites.chr22.20M_filt_gnomADg_AF_anno.vcf.gz  -c INFO/gnomADg_AF -Oz -o gnomad.exomes.v4.1.1.sites.chr22.20M_filt_gnomADg_AF_anno2.vcf.gz  gnomad.exomes.v4.1.1.sites.chr22.20M_filt.vcf.gz
tabix gnomad.exomes.v4.1.1.sites.chr22.20M_filt_gnomADg_AF_anno2.vcf.gz
bcftools query -H  -f '%CHROM\t%POS\t%REF\t%ALT\t%AF\t%gnomADg_AF\n' gnomad.exomes.v4.1.1.sites.chr22.20M_filt_gnomADg_AF_anno2.vcf.gz > testtab.tsv

Then filter for sites present at WGS AF > 0.00002 in R:

chr22 <- read.table('testtab.tsv', header = T, comment.char = '', sep='\t') 
colnames(chr22) <- c("CHROM", "POS","REF","ALT","WES_AF", "WGS_AF")
min_AF <- 0.00002
chr22$WES_AF <- as.numeric(chr22$WES_AF)
chr22$WGS_AF <- as.numeric(chr22$WGS_AF)
chr22_f <- chr22[ !is.na(chr22$WGS_AF) & !is.na(chr22$WES_AF),]
chr22_f <- chr22_f[chr22_f$WES_AF > min_AF & chr22_f$WGS_AF > min_AF,]
hist(log10(chr22_f$WGS_AF / chr22_f$WES_AF) ,1000, xlim=c(-2,2))

Although most AFs do correlate (centering on zero) I can see a sub-population of WGS AFs that is ~10x that of WES and I'm not sure what this means and which one should I report for these variants?

Here are some examples of these variants:

       CHROM      POS REF ALT      WES_AF     WGS_AF
140041 chr22 18855887   C   A 8.38724e-05 0.00054932
35500  chr22 16473030   T   A 1.41497e-04 0.00138468
23493  chr22 15690717   G   T 7.68869e-05 0.00088238
122749 chr22 18370816   G   A 5.32174e-04 0.00462293
98402  chr22 17817297   G   A 2.43971e-04 0.00207800