SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
Understanding the qmap bis-seq format gwilson Epigenetics 3 10-31-2016 12:47 PM
VCF to Circos format(s) bzdyelnik Genomic Resequencing 4 05-27-2014 08:23 AM
problems understanding pileup format pi101 Bioinformatics 2 11-14-2012 02:47 PM
Understanding BAM format. Joker!sAce Genomic Resequencing 7 03-16-2011 06:55 PM

Reply
 
Thread Tools
Old 02-09-2011, 03:02 AM   #1
ketan_bnf
Member
 
Location: India

Join Date: Oct 2010
Posts: 59
Smile Understanding VCF format

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,

Quote:
##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?

Last edited by ketan_bnf; 02-09-2011 at 03:07 AM.
ketan_bnf is offline   Reply With Quote
Old 02-09-2011, 10:27 AM   #2
swbarnes2
Senior Member
 
Location: San Diego

Join Date: May 2008
Posts: 912
Default

"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.
swbarnes2 is offline   Reply With Quote
Old 02-09-2011, 07:58 PM   #3
ketan_bnf
Member
 
Location: India

Join Date: Oct 2010
Posts: 59
Default

yes i have read samtools page.

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

thanks.
ketan_bnf is offline   Reply With Quote
Old 02-09-2011, 11:54 PM   #4
Bruins
Member
 
Location: Groningen

Join Date: Feb 2010
Posts: 78
Default

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.
Bruins is offline   Reply With Quote
Old 02-10-2011, 05:16 AM   #5
laura
Senior Member
 
Location: Cambridge UK

Join Date: Sep 2008
Posts: 151
Default

Quote:
Originally Posted by ketan_bnf View Post
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.
laura is offline   Reply With Quote
Old 02-11-2011, 10:26 AM   #6
Michael.James.Clark
Senior Member
 
Location: Palo Alto

Join Date: Apr 2009
Posts: 213
Default

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.
__________________
Mendelian Disorder: A blogshare of random useful information for general public consumption. [Blog]
Breakway: A Program to Identify Structural Variations in Genomic Data [Website] [Forum Post]
Projects: U87MG whole genome sequence [Website] [Paper]
Michael.James.Clark is offline   Reply With Quote
Old 02-11-2011, 11:41 PM   #7
laura
Senior Member
 
Location: Cambridge UK

Join Date: Sep 2008
Posts: 151
Default

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
laura is offline   Reply With Quote
Old 02-12-2011, 01:00 AM   #8
ketan_bnf
Member
 
Location: India

Join Date: Oct 2010
Posts: 59
Default

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.
ketan_bnf is offline   Reply With Quote
Old 02-12-2011, 01:53 AM   #9
laura
Senior Member
 
Location: Cambridge UK

Join Date: Sep 2008
Posts: 151
Default

Yep bcftools is pretty certain that that individual is homozygous non ref at that position
laura is offline   Reply With Quote
Old 02-21-2011, 03:51 PM   #10
giverny
Member
 
Location: Canada

Join Date: Mar 2010
Posts: 10
Default

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
giverny is offline   Reply With Quote
Old 02-21-2011, 07:16 PM   #11
ketan_bnf
Member
 
Location: India

Join Date: Oct 2010
Posts: 59
Default

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

Last edited by ketan_bnf; 02-21-2011 at 08:30 PM.
ketan_bnf is offline   Reply With Quote
Old 02-22-2011, 07:55 AM   #12
giverny
Member
 
Location: Canada

Join Date: Mar 2010
Posts: 10
Default

Thanks a lot
giverny is offline   Reply With Quote
Old 02-22-2011, 11:25 PM   #13
Bruins
Member
 
Location: Groningen

Join Date: Feb 2010
Posts: 78
Default

Thanks for all the above posts, very helpful!
Bruins is offline   Reply With Quote
Old 04-14-2011, 01:56 AM   #14
marcela
Junior Member
 
Location: Sweden

Join Date: Feb 2011
Posts: 7
Question

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!

Quote:
Originally Posted by ketan_bnf View Post
Hi giverny

PL:GT:GQ 96,47,70,35,0,70:1/1:72
...
P(D|AA) = 10^(-0.96) = 0.109
marcela is offline   Reply With Quote
Old 04-14-2011, 08:02 PM   #15
ketan_bnf
Member
 
Location: India

Join Date: Oct 2010
Posts: 59
Default

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.
ketan_bnf is offline   Reply With Quote
Old 04-19-2011, 06:56 PM   #16
rururara
Member
 
Location: montreal

Join Date: Jan 2011
Posts: 31
Default

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?
rururara is offline   Reply With Quote
Old 04-22-2011, 01:03 PM   #17
Michael.James.Clark
Senior Member
 
Location: Palo Alto

Join Date: Apr 2009
Posts: 213
Default

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)."
__________________
Mendelian Disorder: A blogshare of random useful information for general public consumption. [Blog]
Breakway: A Program to Identify Structural Variations in Genomic Data [Website] [Forum Post]
Projects: U87MG whole genome sequence [Website] [Paper]
Michael.James.Clark is offline   Reply With Quote
Old 05-03-2011, 10:54 AM   #18
blackgore
Member
 
Location: UK

Join Date: Sep 2009
Posts: 20
Default

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.


Quote:
##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?)?

Quote:
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]
blackgore is offline   Reply With Quote
Old 05-03-2011, 08:58 PM   #19
ketan_bnf
Member
 
Location: India

Join Date: Oct 2010
Posts: 59
Default

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.
ketan_bnf is offline   Reply With Quote
Old 05-04-2011, 01:32 AM   #20
blackgore
Member
 
Location: UK

Join Date: Sep 2009
Posts: 20
Default

Hi ketan,

thanks for your fast response! Unfortunately the AF1 tag is not very informative to me, as it's only for the first ALT allele, and pretty much always 1.0 or 0.5. The AF tag, as I understand it, gives proportional representation of all the called alleles, which is what I'm really after.
blackgore is offline   Reply With Quote
Reply

Tags
vcf

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off




All times are GMT -8. The time now is 12:15 AM.


Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2019, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO