[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?

[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

[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

[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!

[aforntacc]

OK, i have opened a new thread,
i have only one objective to compare wt (0/1:93,0,255:33:0:95, 0/0:0,122,255:59:0:99 0/1:53,0,241:42:0:55) and mutant (0/1:139,0,250:59:7:99, 0/1:78,0,255:70:3:80 0/1:62,0,238:28:0:64).
to find what is common indels and unique ones
and i kindly ask please do i interpret the indels as the way the snps are interpreted. as in ref =0 alt =1. because i tried to read stuff about indels i cant find anything useful so i turn to this site for the much need help.
thanks a lot for your help i am grateful

[laura]

You may find it more fruitful to start another post with this question with a more specific title

How do you interpret a particular indel is really a how long is a piece of string question, what are you trying to find out, what hypothesis are you testing and these will determine how you might look at your result

vcf tools has various scripts to produce intersections and compare two different vcf files, which tools are appropriate again very much depends on what questions you want to answer

I would think about the specific questions you wish to answer and come back with a new thread and more information and you may get more help that way

[aforntacc]

Please i want to find out how i can interpret the indels, i have some thing like this chr3 4466963 . TGGAG TGGAGGAG 999 PASS AC1=5;AF1=0.4167;DP4=7,197,3,84;DP=347;FQ=999;G3=0.1667,0.8333,8.319e-50;HWE=0.0465;INDEL;MQ=44;MfGt=0/1;MinDP=28;NeqMfGt=1;PV4=1,1.5e-70,2e-112,1 GT:PL: DP:SP:GQ 0/1:93,0,255:33:0:95 0/0:0,122,255:59:0:99 0/1:53,0,241:42:0:55 0/1:139,0,250:59:7:99 0/1:78,0,255:70:3:80 0/1:62,0,238:28:0:64
so will it be like this TGGAG = 0 and TGGAGGAG = 1 as the same for the snps. or is it interpreted differently?
secondly how can i use the vcftols to compare two individuals (i have seperated them into two vcf files) with 3 libraries WT (0/1:93,0,255:33:0:95 0/0:0,122,255:59:0:99 0/1:53,0,241:42:0:55) and
MT (0/1:139,0,250:59:7:99 0/1:78,0,255:70:3:80 0/1:62,0,238:28:0:64).
thanks a lot[/QUOTE]

[aforntacc]

Hello

Please i want to find out how i can interpret the indels, i have some thing like this chr3 4466963 . TGGAG TGGAGGAG 999 PASS AC1=5;AF1=0.4167;DP4=7,197,3,84;DP=347;FQ=999;G3=0.1667,0.8333,8.319e-50;HWE=0.0465;INDEL;MQ=44;MfGt=0/1;MinDP=28;NeqMfGt=1;PV4=1,1.5e-70,2e-112,1 GT:PL: DP:SP:GQ 0/1:93,0,255:33:0:95 0/0:0,122,255:59:0:99 0/1:53,0,241:42:0:55 0/1:139,0,250:59:7:99 0/1:78,0,255:70:3:80 0/1:62,0,238:28:0:64
so will it be like this TGGAG = 0 and TGGAGGAG = 1 as the same for the snps. or is it interpreted differently?
secondly how can i use the vcftols to compare two individuals (i have seperated them into two vcf files) with 3 libraries WT (0/1:93,0,255:33:0:95 0/0:0,122,255:59:0:99 0/1:53,0,241:42:0:55) and
MT (0/1:139,0,250:59:7:99 0/1:78,0,255:70:3:80 0/1:62,0,238:28:0:64).
thanks a lot

[laura]

The missing ids are unlikely to be the cause of this

You need to provide specific details about the filtering methods you are choosing if people are going to help figure out what is wrong

[aforntacc]

Hi all

i dont have snp id in my data, the ID column is all in dot (.) why is this because i am able to filter out the indels but not the snps
how can i do this.
thanks
#CHROM POS ID REF ALT
chr1 8686 . T C
chr1 10802 . T C
chr1 10815 . A G
chr1 10836 . C A
chr1 11355 . C A
chr1 11433 . T C
chr1 11669 . ATTTT ATTTTT

[laura]

It looks like you either need to investigate if your problem can be solved with the vcftools perl scripts or maybe change your header from vcf4.1 to vcf4.0 and see what the vcftools binary does

These questions are now most appropriate for the vcftools-help list which you can find

http://vcftools.sourceforge.net/

[aforntacc]

hi laura
i am still progressing small small
but i have got this error when i want to output the vcf file with passed snps ( --remove-filtered-geno-all)
bilbo@ubuntu:~/vcftools_0.1.4a$ ./cpp/vcftools --vcf /media/My\ Passport/other\ analysis\ by\ fasteris/2012-02-21_GQJ-1-6_VitisVinifera_variants.vcf --remove-filtered-geno-all --out /media/My\ Passport/other\ analysis\ by\ fasteris/lagolas.vcf

VCFtools - v0.1.4
(C) Adam Auton 2009

Parameters as interpreted:
--out /media/My Passport/other analysis by fasteris/lagolas.vcf
--remove-filtered-geno-all
--vcf /media/My Passport/other analysis by fasteris/2012-02-21_GQJ-1-6_VitisVinifera_variants.vcf

Scanning /media/My Passport/other analysis by fasteris/2012-02-21_GQJ-1-6_VitisVinifera_variants.vcf ...
Error:VCF version must be v4.0:
You are using version VCFv4.1

now i am stuck, please what should i do.
thanks

[laura]

Unfortunately that is a bit of a how long is a piece of string question as it very much depends on what tools/programming language you wish to use to do it

If you want a vcf file with just PASS snps in it you can use the vcftools binary and its --remove-filtered-geno-all option but if you want other info than that then it depends

[aforntacc]

Thank you Laura, gradually i am making progress.
please i want to ask but i dont know if this is a stupid question
if i want to uncode the GT index for all SNPs that pass the filter criteria how can i do that?
specifically do i have to do this with the VCF tools (decode genotype) using the PERL5LIB environment or what? am a bit confused please.
thanks a lot

[laura]

If you want to know the difference between your WT and MT individuals you need to compare their genotypes