How to Recalculate Allele Frequency from a VCF File?

Question

I have a VCF file with 200 samples (mitochondrial genome of Plasmodium falciparum). Here is a pic to take a look at:

And a few relevant lines from the actual file:

##INFO=<ID=AC,Number=A,Type=Integer,Description="Allele count in genotypes, for each ALT allele, in the same order as listed">
##INFO=<ID=AF,Number=A,Type=Float,Description="Allele Frequency, for each ALT allele, in the same order as listed">
##INFO=<ID=AN,Number=1,Type=Integer,Description="Total number of alleles in called genotypes">
#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  FP0008-C   FP0009-C
Pf_M76611       19      .       A       T       33196.8 MissingVQSLOD;Mitochondrion     AC=9;AF=0.0003019;AN=29816;ANN=T|intragenic_variant|MODIFIER|PF3D7_MIT00100|PF3D7_MIT00100|gene_variant|PF3D7_MIT00100|||n.19A>T||||||;AS_InbreedingCoeff=0.0437;AS_QD=20.96;BaseQRankSum=0.141;CDS;DP=559138;ExcessHet=0.0026;FS=2.892;InbreedingCoeff=0.52;MLEAC=8;MLEAF=0.0002683;MQ=60;MQRankSum=0;QD=20.96;ReadPosRankSum=0.856;RegionType=Mitochondrion;SOR=0.789;set=genotypegvcfs  GT:AD:DP:GQ:PL  0/0:52,0:52:99:0,120,1800  0/0:54,0:54:99:0,120,1800

A small subset of samples I am working on can be found here: https://www.mediafire.com/file/6dkct6guk5zx83m/sample.vcf/file

I am trying to calculate allele frequency (AF) for each variant in the dataset. The population is assumed to be all the samples in the file.

In other words, even though the VCF file contains AF statistics from population data, I would like to know how to recalculate AF only from the VCF file. I believe I need to know how AC and AN are calculated first.

AF is calculated this way: allele count (AC) / allele number (AN)

I tried to calculate them from the information from GT (counting 0/0, 1/1 and 0/1 in all the 200 samples), but then I learned that this calculation is not suitable for mithocondrial data.

All lines in the VCF file contain an AF statistic:

$ grep -v 'AF=' sample.vcf | grep -v '^#'
<No output>

The information is provided under INFO column, but how were they calculated? If I do not have AC and AN to begin with, what should I do? Is there a way to calculate AC and AN from the other information provided in the dataset? such as the data mentioned in the FORMAT column? GT: genotype AD: unfiltered allele depth DP : read depth at this position for this sample (Integer) GQ : conditional genotype quality, encoded as a phred quality PL : the phred-scaled genotype likelihoods rounded to the closest integer

Thanks

Answer 1

It sounds like your main confusion is how to calculate AF. I'll answer that right away, but I'll also suggest a Python library that you can use to do this and much more.

You're correct, the AF is defined as the number of non-missing alternate alleles (AC) divided the number of non-missing alleles (AN). Calculating these are straightforward. Map each GT to its alternate alleles:

0/0 => 0
0/1 => 1
1/1 => 2
./0 => 0
./1 => 1
./. => 0

Now, compute the sum. That is the AC. For the AN, map each GT to its number of non-missing alleles:

0/0 => 2
0/1 => 2
1/1 => 2
./0 => 1
./1 => 1
./. => 0

Now, compute the sum. That is the AN. Now calculate AF as

AC / AN

All of this is kind of annoying to do though! I suggest you use Hail (disclaimer: I am a Hail developer).

import hail as hl

# Hail is not yet aware of Plasmodium falciparum reference genome
mt = hl.import_vcf('sample.vcf', reference_genome=None)

# let's peek at our data
mt.show(n_rows=10, n_cols=5)

# Let's get some information about the GT specifically:
mt.GT.summarize()

# We can use Hail's built in variant_qc method to calculate AF for us:
mt = hl.variant_qc(mt)

# Let's take a look at that
mt.variant_qc.show()

# We can overwrite the old AF field and export back to a VCF file if necessary
mt = mt.annotate_rows(
    # AF[1] is the frequency of the first alternate allele
    # Your dataset appears to have some multi-allelics. The other
    # alternate allele frequencies are available in AF[2], ...
    info = mt.info.annotate(
        AF = mt.variant_qc.AF[1]
    )
)
hl.export_vcf(mt, 'sample-updated.vcf.bgz')

---

This comment doesn't make sense to me:

I learned that this calculation is not suitable for mithocondrial data.

The definition of AF, AC, and AN doesn't change depending on the type of config (autosome, sex chromosome, mitochondrial). Perhaps the issue is ploidy encoding? A GT in a VCF can be:

• 1-ploidy, e.g. 0 (reference), 1 (first alternate), 2 (second alternate)
• 2-ploidy, e.g. 0/0 (homozygous reference), 0/1 (heterozygous), 1/1 (homozygous first alternate)
• etc.

Sometimes people incorrectly encode a 1-ploidy config (for example the Y chromosome) as: 0/0 and 1/1 instead of 0 and 1. That's not your fault though and it doesn't mean the definition of AF, AC, and AN change. It means someone gave you malformed data.