[ketan_bnf]

Hi blackgore,

for your first question about Allele frequency, in vcf file AF1 is showing the Allele frequency for the related alt base.

read starting of the vcf file there are coments,

##INFO=<ID=AF1,Number=1,Type=Float,Description="Max-likelihood estimate of the site allele frequency of the first ALT allele">

chr2.fa 102741718 . A G,C,T 14.2 . DP=5511;AF1=0.5;CI95=0.5,0.5;DP4=2024,2411,432,528;MQ=58;FQ=17.1;PV4=0.75,3.5e-06,0.0085,1 GT:PL:GQ 0/1:44,0,248,255,255,255,255,255,255,255:47

for your 2nd question i am not sure.

[blackgore]

Hi,

At the risk of appearing quite dumb, if anyone can help with (admittedly my first day using) the vcf format and vcftools, I'd be very grateful! For reference, I'm using illumina reads, samtools v0.1.13, and vcftools 0.1.15. BAM alignment was created using BWA.

VCF file was generated by the following:

samtools view -b [bamFile] [regions of interest] | mpileup -uf [reference genome] - | bcftools view -vcgAN - > variants.raw.bcf

and validated using

vcf-validator variants.raw.bcf

My questions (so far) are:

1) How I can get AF (allele frequency) values into the VCF file? Despite how I've tried, it does not seem to want to appear in the output.

##fileformat=VCFv4.1
##samtoolsVersion=0.1.13 (r926:134)
##INFO=<ID=DP,Number=1,Type=Integer,Description="Raw read depth">
##INFO=<ID=DP4,Number=4,Type=Integer,Description="# high-quality ref-forward bases, ref-reverse, alt-forward and alt-reverse bases">
##INFO=<ID=MQ,Number=1,Type=Integer,Description="Root-mean-square mapping quality of covering reads">
##INFO=<ID=FQ,Number=1,Type=Float,Description="Phred probability that sample chromosomes are not all the same">
##INFO=<ID=AF1,Number=1,Type=Float,Description="Max-likelihood estimate of the site allele frequency of the first ALT allele">
##INFO=<ID=CI95,Number=2,Type=Float,Description="Equal-tail Bayesian credible interval of the site allele frequency at the 95% level">
##INFO=<ID=PV4,Number=4,Type=Float,Description="P-values for strand bias, baseQ bias, mapQ bias and tail distance bias">
##INFO=<ID=INDEL,Number=0,Type=Flag,Description="Indicates that the variant is an INDEL.">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##FORMAT=<ID=GQ,Number=1,Type=Integer,Description="Genotype Quality">
##FORMAT=<ID=GL,Number=3,Type=Float,Description="Likelihoods for RR,RA,AA genotypes (R=ref,A=alt)">
##FORMAT=<ID=DP,Number=1,Type=Integer,Description="# high-quality bases">
##FORMAT=<ID=SP,Number=1,Type=Integer,Description="Phred-scaled strand bias P-value">
##FORMAT=<ID=PL,Number=-1,Type=Integer,Description="List of Phred-scaled genotype likelihoods, number of values is (#ALT+1)*(#ALT+2)/2">
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT -
chr2.fa 102741718 . A G,C,T 14.2 . DP=5511;AF1=0.5;CI95=0.5,0.5;DP4=2024,2411,432,528;MQ=58;FQ=17.1;PV4=0.75,3.5e-06,0.0085,1 GT:PL:GQ 0/1:44,0,248,255,255,255,255,255,255,255:47
chr2.fa 102742168 . T C,G,A 143 . DP=5872;AF1=0.5;CI95=0.5,0.5;DP4=1633,2355,778,917;MQ=59;FQ=146;PV4=0.0006,1,2.2e-21,1 GT:PL:GQ 0/1:173,0,255,255,255,255,255,255,255,255:99
chr2.fa 102745722 . AAC AACNAC 217 . INDEL;DP=2851;AF1=0.5;CI95=0.5,0.5;DP4=357,8,604,10;MQ=59;FQ=217;PV4=0.62,1,0.0028,0.12 GT:PL:GQ 0/1:255,0,255:99
chr2.fa 102748084 . G A,X 5.46 . DP=22;AF1=0.4999;CI95=0.5,0.5;DP4=10,6,2,3;MQ=59;FQ=7.8;PV4=0.61,0.00013,0.036,1 GT:PL:GQ 0/1:34,0,255,82,255,255:37
chr2.fa 102748094 . A G,X 21 . DP=23;AF1=0.5;CI95=0.5,0.5;DP4=13,5,2,3;MQ=59;FQ=24;PV4=0.3,0.00017,0.028,1 GT:PL:GQ 0/1:51,0,255,105,255,255:54
chr2.fa 102750361 . CAAAAAAAAAAAAA CAAAAAAAAAAA,CAAAAAAAAAAAAAA 29 . INDEL;DP=38;AF1=1;CI95=0.5,1;DP4=10,0,20,4;MQ=54;FQ=-46.5;PV4=0.3,1,1,1 GT:PL:GQ 1/1:69,12,0,83,39,73:72
chr2.fa 102750722 . GTG GTGTTCTCTG,GTTCTCTG 217 . INDEL;DP=3246;AF1=0.5;CI95=0.5,0.5;DP4=374,410,913,997;MQ=42;FQ=217;PV4=0.97,1,0,1 GT:PL:GQ 0/1:255,0,255,255,255,255:99
chr2.fa 102751604 . GTTTTTTTT GTTTTTTT,GTTTTTTTTT 56.5 . INDEL;DP=5535;AF1=1;CI95=1,1;DP4=456,504,1961,2049;MQ=59;FQ=-246;PV4=0.45,1,5.9e-10,0.24 GT:PL:GQ 1/1:97,211,0,244,255,123:99

2) The validator script's output (snippet below) gives me a number of issues, which I'd like to sort out. How do I add the missing header information? Without knowing enough about the vcf format just yet, I'm shying away from manually opening the file and adding the information. Also, does the 'X' allele signify a problem, or is it just an "unknown" allele (but why would it be?)?

The header tag 'reference' not present. (Not required but highly recommended.)
The header tag 'contig' not present for CHROM=chr2.fa. (Not required but highly recommended.)
chr2.fa:102740231 .. Could not parse the allele(s) [X]
chr2.fa:102740459 .. Could not parse the allele(s) [X]
chr2.fa:102748084 .. Could not parse the allele(s) [X]
chr2.fa:102748094 .. Could not parse the allele(s) [X]

[Michael.James.Clark]

Hi rururara

Read the VCF4 specs here:

http://www.1000genomes.org/node/101

"If genotype information is present, then the same types of data must be present for all samples. First a FORMAT field is given specifying the data types and order. This is followed by one field per sample, with the colon-separated data in this field corresponding to the types specified in the format. The first sub-field must always be the genotype (GT)."

[rururara]

hai all,

i try to validate my snp in vcf format using vcf-validator.
But it print out as below :


Expected GT as the first genotype field at Chr12:15841735
Expected GT as the first genotype field at Chr12:17651041
Expected GT as the first genotype field at Chr12:17804331
Expected GT as the first genotype field at Chr12:18935754
Expected GT as the first genotype field at Chr12:19270259
Expected GT as the first genotype field at Chr12:19878395
Expected GT as the first genotype field at Chr12:19951137


Can somebody explain me what does it means? I try to find the information but still not clear.
Is that an error or anything?

Thanks?

[ketan_bnf]

Hi! marcela,

Thanks for pointing out an error, my mistake in writing PL,

it should be like

(...) A T,G (...) DP=35;AF1=1;CI95=0.5,1;DP4=0,0,1,8;MQ=60 PL:GT:GQ 96,47,70,35,0,70:1/1:72

P(D|AA) = 10^(-9.6) = ?

like on wards,

Thanks.

[marcela]

Hi!

just a question, at http://samtools.sourceforge.net/mpileup.shtml, there is this example:

REF=C
ALT=A,G,

PL=7,0,37,13,40,49

P(D|CC)=10^{-0.7}
P(D|CA)=1
P(D|AA)=10^{-3.7}
P(D|CG)=10^{-1.3}
P(D|AG)=1e-4
P(D|GG)=10^{-4.9}

From yours, when calculating P(D|AA) you use 0.96 instead of 9.6 if we follow the example above.. which is correct?

Thanks!

Hi giverny

PL:GT:GQ 96,47,70,35,0,70:1/1:72
...
P(D|AA) = 10^(-0.96) = 0.109

[Bruins]

Thanks for all the above posts, very helpful!

[giverny]

Thanks a lot

[ketan_bnf]

Hi giverny

(...) A T,G (...) DP=35;AF1=1;CI95=0.5,1;DP4=0,0,1,8;MQ=60 PL:GT:GQ 96,47,70,35,0,70:1/1:72

ref = A, alt = T,G right?

So, genotype may be AA,AT,AG,TT,TG,GG

PL means phread-scaled likelihood

P(D|AA) = 10^(-0.96) = 0.109
P(D|AT) = 10^(-0.47) = 0.33
P(D|AG) = 10^(-0.70) = 0.199
P(D|TT) = 10^(-0.35) = 0.446
P(D|TG) = 10^(0) = 1
P(D|GG) = 10^(-0.70) = 0.199

So, according to your genotype likelihood, genotype is TG, if the genotype what i have count is right.

you can visit http://samtools.sourceforge.net/mpileup.shtml , SAMtools/BCFtools specific information

[giverny]

Hi there,
I have more than one alternative allele (see below)

(...) A T,G (...) DP=35;AF1=1;CI95=0.5,1;DP4=0,0,1,8;MQ=60 PL:GT:GQ 96,47,70,35,0,70:1/1:72

In this case, how do I have to read the likelihood please? Thanks a lot

[laura]

Yep bcftools is pretty certain that that individual is homozygous non ref at that position

[ketan_bnf]

chr1 10740313 . A G 188.30 PASS AC=2;AF=1.00;AN=2;DP=11;Dels=0.00;HRun=1;Haplotype Score=6.9635;MQ=26.82;MQ0=0;QD=17.12;SB=-72.04;sumGLbyD=20.12 GT:AD: DP :GQ:PL 1/1:1,10:7:21.05:221,21,0

Here PL is 221,21,0

according to samtools mpileup page

PL means SAMtools/BCFtools writes genotype likelihoods in the PL format which is a comma delimited list of phred-scaled data likelihoods of each possible genotype.

P(D|AA) = 10^(-2.21) = 0.006
P(D|AG) = 10^(-0.21) = 0.617
P(D|GG) = 10^(0) = 1

so does it means genotype is GG for this SNP?

And thanks for AD and DP, now i understood it.

[laura]

The important thing to under stand any VCF file is the read the header and if that isn't clear enough go to the documentation from the caller. Why the GATK developers have chosen to represent AD and total DP differently I am not sure but I do know if you go and ask here http://getsatisfaction.com/gsa they will probably be able to explain

[Michael.James.Clark]

Something that has thrown myself and others (I saw this same issue on the Biostar Stackexchange) for a loop is that "sum(AD) ne DP". It would seem that the sum of the allelic depths should equal the total depth, but this is not the case from GATK at least--it's always equal to or greater than the total depth.

The reason seems to be because the DP is only filtered reads (reads used for calling) whereas AD is (presumably) all reads.

Of course, I'm not sure why we'd keep unfiltered info that wasn't used for calling in the first place there, though I think there must be a rationale. Still, it makes me question using AD to assess reference calling bias generally.

[laura]

Hi all,

I need help to understand vcf format. i have visited all sites (1000genomes, vcf specifications etc) still need some understanding,

Called SNPs with GATK,



GT means genotype but what is phased and unphased? what does 0/1 or 1/1 means for GT?
what is the meaning of PL, HaplotypeScore, SB?

Is AF depends on AC and AN value?

Should DP (11) of INFO tag and sum of DP (1,10:ref reads, alt reads) of FORMAT tag always be the same?

First have you read the documentation about the format which you can find here

http://vcftools.sourceforge.net/specs.html

The header of the vcf file gives you some information about each field, if you were using GATK the GATK documentation should also explain a bit

As far as your specific questions go AF will depend on AC and AN values but if LD information was used to calculate the AF it won't necessarily be the same as AC/AN

For Depth I suspect it depends on how the caller calculated depth for each individual and total depth but I would be surprised if DP in the info column wasn't the sum of DP in the individual fields because in both instances its mean to represent the depth of reads used to call a particular variant.