Understanding VCF format

[ketan_bnf]

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,

##fileformat=VCFv4.0
##FILTER=<ID=LowQual,Description="Low quality">
##FORMAT=<ID=AD,Number=.,Type=Integer,Description="Allelic depths for the ref and alt alleles in the order listed">
##FORMAT=<ID=DP,Number=1,Type=Integer,Description="Read Depth (only filtered reads used for calling)">
##FORMAT=<ID=GQ,Number=1,Type=Float,Description="Genotype Quality">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##FORMAT=<ID=PL,Number=3,Type=Float,Description="Normalized, Phred-scaled likelihoods for AA,AB,BB genotypes where A=ref and B=alt; not applicable if site is not biallelic">
##INFO=<ID=AC,Number=.,Type=Integer,Description="Allele count in genotypes, for each ALT allele, in the same order as listed">
##INFO=<ID=AF,Number=.,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">
##INFO=<ID=DP,Number=1,Type=Integer,Description="Total Depth">
##INFO=<ID=DS,Number=0,Type=Flag,Description="Were any of the samples downsampled?">
##INFO=<ID=Dels,Number=1,Type=Float,Description="Fraction of Reads Containing Spanning Deletions">
##INFO=<ID=HRun,Number=1,Type=Integer,Description="Largest Contiguous Homopolymer Run of Variant Allele In Either Direction">
##INFO=<ID=HaplotypeScore,Number=1,Type=Float,Description="Consistency of the site with two (and only two) segregating haplotypes">
##INFO=<ID=MQ,Number=1,Type=Float,Description="RMS Mapping Quality">
##INFO=<ID=MQ0,Number=1,Type=Integer,Description="Total Mapping Quality Zero Reads">
##INFO=<ID=QD,Number=1,Type=Float,Description="Variant Confidence/Quality by Depth">
##INFO=<ID=SB,Number=1,Type=Float,Description="Strand Bias">
##UnifiedGenotyper="analysis_type=UnifiedGenotyper input_file=[../../samtools-0.1.12a/samfiles/q20/GKUNU9Q04_chr1_q20_RG_sort.bam] sample_metadata=[] read_buffer_size=null phone_home=STANDARD read_filter=[] intervals=null excludeIntervals=null reference_sequence=chr1.fa rodBind=[] rodToIntervalTrackName=null BTI_merge_rule=UNION DBSNP=null downsampling_type=null downsample_to_fraction=null downsample_to_coverage=50 baq=OFF baqGapOpenPenalty=40.0 performanceLog=null useOriginalQualities=false defaultBaseQualities=-1 validation_strictness=SILENT unsafe=null num_threads=1 interval_merging=ALL read_group_black_list=null processingTracker=null restartProcessingTracker=false processingTrackerStatusFile=null processingTrackerID=-1 logging_level=INFO log_to_file=null quiet_output_mode=false debug_mode=false help=false genotype_likelihoods_model=SNP p_nonref_model=EXACT heterozygosity=0.001 pcr_error_rate=1.0E-4 sites_only=false genotype=false output_all_callable_bases=false standard_min_confidence_threshold_for_calling=50.0 standard_min_confidence_threshold_for_emitting=10.0 trigger_min_confidence_threshold_for_calling=30.0 trigger_min_confidence_threshold_for_emitting=30.0 noSLOD=false assume_single_sample_reads=null min_base_quality_score=17 min_mapping_quality_score=20 max_mismatches_in_40bp_window=3 use_reads_with_bad_mates=false max_deletion_fraction=0.05 get_indel_alleles_from_vcf=false min_indel_count_for_genotyping=5 indel_heterozygosity=1.25E-4 out=org.broadinstitute.sting.gatk.io.stubs.VCFWriterStub NO_HEADER=org.broadinstitute.sting.gatk.io.stubs.VCFWriterStub verbose_mode=null metrics_file=null annotation=[]"
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT buffalo

chr1 10740313 . A G 188.30 PASS AC=2;AF=1.00;AN=2;DP=11;Dels=0.00;HRun=1;HaplotypeScore=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

chr1 147668098 . A T 12.05 LowQual AC=1;AF=0.50;AN=2;DP=6;Dels=0.00;HRun=0;HaplotypeScore=3.7429;MQ=121.66;MQ0=0;QD=2.01;SB=-0.01;sumGLbyD=7.25 GT:AD: DP: GQ :PL 0/1:0,6:1:1.76:42,3,0

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?

[swbarnes2]

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

Not necessarily. Here's the quote from the samtools page:

Number of 1) forward ref alleles; 2) reverse ref; 3) forward non-ref; 4) reverse non-ref alleles, used in variant calling. Sum can be smaller than DP because low-quality bases are not counted.

[ketan_bnf]

yes i have read samtools page.

if anybody can give more insights that will be more helpfull.

thanks.

[Bruins]

I struggle(d) with this too. Googled much, found little answers. Perhaps an idea to add this kind of thing (where to find docs on file formats + additional explanation) to the wiki? It would help to have a central point to start from.

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

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

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

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

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

[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

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

Thanks a lot

[Bruins]

Thanks for all the above posts, very helpful!

[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

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