Seqanswers Leaderboard Ad

Collapse

Announcement

Collapse
No announcement yet.
X
 
  • Filter
  • Time
  • Show
Clear All
new posts

  • 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,

    ##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, 04:07 AM.

  • #2
    "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.

    Comment


    • #3
      yes i have read samtools page.

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

      thanks.

      Comment


      • #4
        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.

        Comment


        • #5
          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



          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.

          Comment


          • #6
            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]

            Comment


            • #7
              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

              Comment


              • #8
                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.

                Comment


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

                  Comment


                  • #10
                    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

                    Comment


                    • #11
                      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, 09:30 PM.

                      Comment


                      • #12
                        Thanks a lot

                        Comment


                        • #13
                          Thanks for all the above posts, very helpful!

                          Comment


                          • #14
                            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!

                            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

                            Comment


                            • #15
                              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.

                              Comment

                              Latest Articles

                              Collapse

                              • seqadmin
                                Strategies for Sequencing Challenging Samples
                                by seqadmin


                                Despite advancements in sequencing platforms and related sample preparation technologies, certain sample types continue to present significant challenges that can compromise sequencing results. Pedro Echave, Senior Manager of the Global Business Segment at Revvity, explained that the success of a sequencing experiment ultimately depends on the amount and integrity of the nucleic acid template (RNA or DNA) obtained from a sample. “The better the quality of the nucleic acid isolated...
                                03-22-2024, 06:39 AM
                              • seqadmin
                                Techniques and Challenges in Conservation Genomics
                                by seqadmin



                                The field of conservation genomics centers on applying genomics technologies in support of conservation efforts and the preservation of biodiversity. This article features interviews with two researchers who showcase their innovative work and highlight the current state and future of conservation genomics.

                                Avian Conservation
                                Matthew DeSaix, a recent doctoral graduate from Kristen Ruegg’s lab at The University of Colorado, shared that most of his research...
                                03-08-2024, 10:41 AM

                              ad_right_rmr

                              Collapse

                              News

                              Collapse

                              Topics Statistics Last Post
                              Started by seqadmin, Yesterday, 06:37 PM
                              0 responses
                              11 views
                              0 likes
                              Last Post seqadmin  
                              Started by seqadmin, Yesterday, 06:07 PM
                              0 responses
                              10 views
                              0 likes
                              Last Post seqadmin  
                              Started by seqadmin, 03-22-2024, 10:03 AM
                              0 responses
                              51 views
                              0 likes
                              Last Post seqadmin  
                              Started by seqadmin, 03-21-2024, 07:32 AM
                              0 responses
                              68 views
                              0 likes
                              Last Post seqadmin  
                              Working...
                              X