2
$\begingroup$

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

enter image description here

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

$\endgroup$
2
  • $\begingroup$ Can you give us some more context? What is AC? Allele count in the population? From GnomAD, maybe? What tool produced this file? $\endgroup$
    – terdon
    Jan 27, 2023 at 16:58
  • $\begingroup$ @terdon Thanks for the reply. I believe that is what AC means. I have not seen any other definition. $\endgroup$
    – eh329
    Jan 27, 2023 at 17:56

1 Answer 1

2
$\begingroup$

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.

$\endgroup$
2
  • 1
    $\begingroup$ Thanks for the reply. You totally understood my question. Reagrding the part you said that does not make sense, what i meant was that cells can have multiple mts. From the info i gathered so far, due to this reason the simple equation of AF = AC / AN would not be appropriate for mt genome meaning using GT alone will not return the right result. The calculations used in softwares such as vcftools are only for non mt genome. I was trying to do it in Python anyway. I am glad you introduced this library. Does it work on Windows too? $\endgroup$
    – eh329
    Jan 30, 2023 at 17:04
  • $\begingroup$ We don't regularly test against windows, but, in principle, it should work. Let me know if you run into bugs and we'll fix them! $\endgroup$ Jan 30, 2023 at 21:09

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.