[daisieh]

I think this is relevant to this thread, which is why I'm reawakening it:

According to 1000genomes's VCF 4.1 spec, the ordering of genotypes is given by this:

If A is the allele in REF and B,C,... are the alleles as ordered in ALT, the ordering of genotypes for the likelihoods is given by: F(j/k) = (k*(k+1)/2)+j. In other words, for biallelic sites the ordering is: AA,AB,BB; for triallelic sites the ordering is: AA,AB,BB,AC,BC,CC, etc.

Just in case anyone else is desperately googling for the answer to how to order genotypes for bi/triallelic alternate alleles in a vcf file!

[vyellapa]

Would anyone know why Samtools does not give GT and GQ for some calls. Example:3 52720080 in vcf snippet below.

3 33434831 . G A 14.2 . DP=135;VDB=1.109343e-02;RPB=1.478724e+00;AF1=0.5;AC1=1;DP4=49,47,18,19;MQ=20;FQ=17.1;PV4=0.85,1,0.065,1 GT:PL:GQ 0/1:44,0,113:47
3 42251263 . C T 133 . DP=21;VDB=1.850543e-01;RPB=-1.073440e+00;AF1=1;AC1=2;DP4=0,1,12,8;MQ=20;FQ=-70;PV4=0.43,0.27,1,1 GT:PL:GQ 1/1:166,43,0:83
3 42787469 . A G 9.31 . DP=2;VDB=6.720000e-02;AF1=1;AC1=2;DP4=0,0,1,1;MQ=20;FQ=-33 GT:PL:GQ 1/1:40,6,0:8
3 52720080 . A . 48.9 . DP=67;VDB=1.892600e-02;RPB=-4.838016e-01;AF1=0;AC1=0;DP4=52,0,13,0;MQ=20;FQ=-46;PV4=1,1,1,0.34 PL 0
3 101576175 . T C 171 . DP=100;VDB=1.334901e-01;AF1=1;AC1=2;DP4=0,0,59,40;MQ=20;FQ=-282 GT:PL:GQ 1/1:204,255,0:99

[laura]

vyellapa

This is relevant to samtools and not to the VCF spec. I would recommend creating a new question about this or emailing the samtools help list you will get a better response this way

[drdna]

There are a number of quirks of VCF files that I do not understand. I would appreciate it if someone could explain them to me:

1) VCF reports multiple alternate alleles but provides detailed information for only one of them. As an example, view the following three lines from a samtools mpileup run:

supercont8.1 889760 . T G,C 222 . DP=17;VDB=1.071396e-01;AF1=1;AC1=2;DP4=0,0,13,4;MQ=39;FQ=-75 GT:PL:GQ 1/1:255,48,0,255,37,252:93
supercont8.1 893978 . C T,G 174 . DP=37;VDB=2.882637e-01;AF1=1;AC1=2;DP4=0,0,12,24;MQ=21;FQ=-132 GT:PL:GQ 1/1:207,105,0,207,89,204:99
supercont8.1 905324 . T C,G 213 . DP=44;VDB=2.155955e-01;AF1=1;AC1=2;DP4=0,0,20,16;MQ=27;FQ=-132 GT:PL:GQ 1/1:246,105,0,248,92,245:99

First, the MLE for the first alt allele count is 2. However, if we are only counting the FIRST alt allele, this value should always be 1 (you can only have one FIRST alt allele).

If we disregard this semantic error and assume that there are two alt alleles, then why is the genotype listed as 1/1. After all, in each entry, the data presented thus far suggest two alt alleles. In this case, the genotypes should be 1/2. The list of Phred-scaled genotype likelihoods (PL values) imply low confidence for one the "2" alt allele. Problem, is we are provided no information on the frequency/quality of alt allele #2, so we cannot independently evaluate the mpileup call.

2) Consider the following output line:

supercont8.4 3182995 . AG AACG,AGCCCAACG 30.9 . INDEL;IDV=11;IMF=0.366667;DP=30;VDB=3.034881e-04;AF1=0.8294;AC1=1;DP4=1,0,8,1;MQ=36;FQ=-33.5;PV4=1,1,0.22,0.13 GT:PL:GQ 0/1:120,53,54,67,0,71:4

Here, two alt alleles are reported, however the MLE for first alt allele count is only 1. Presumably mpileup did not like one of the alt allele calls but the vcf report provides no way to find out why.

Is there a way to make mpileup report frequencies/qualities for all alt alleles?

3) Finally, when two alternate alleles are found, I am assuming the list of Phred-scaled genotype likelihoods is provided in the following order: 0/0, 0/1, 0/2, 1/1, 1/2, 2/2. Is this a reasonable assumption?