Unconfigured Ad

Collapse
X
 
  • Filter
  • Time
  • Show
Clear All
new posts
  • jmartin
    Member
    • Dec 2009
    • 78

    Need help interpreting VCF output

    I ran samtools mpileup & then vcfutils varFilter to call snps on a set of 10 illumina samples against a number of reference genomes. I'm now looking at the results, and I am having trouble understanding the information I've been given per sample.

    The VCF file lists my 10 samples as the last 10 columns, and has them in the format GT:PL:GQ. I understand that those show me the genotype, phred-scaled likelihood of the genotypes, and the overall quality of the genotype assignment. But I am confused by the per-sample results. Here are some specific questions I have:

    1) I do not expect every sample to have coverage over every loci called. What should I see in the field for a sample that does not contain that loci? I suspect it looks like this:

    0/1:0,0,0:3

    since that seems to appear a lot in my output. But I want to be sure.

    2) The VCF documentation says that for diploid calls (my samples are actually metagenomic, but the caller seems to consider them diploid so I am going with it), the 'PL' portion of the sample info fields is giving me the phred-scaled likelihood of the genotypes in the order:

    AA AB BB

    Where A is the ref allele & B is the alternate allele. But I am seeing data that looks like this:

    1/1:145,21,0:26

    From the documentation, I had thought a GT value of 1/1 meant that it was seeing, for that given sample, both alleles being the alternate allele. And from the PL documentation, I would have expected the large PL value (145) to be in the THIRD position. I thought the first position in that string was reserved for REF/REF genotype. Could somebody clarify the order of the PL values for me?


    Thanks,
    John Martin

    p.s. The documentation I am referring to is the VCF format v4.1 manual page at the 1000 genomes website (not sure if its OK to give URLs in this forum, but its pretty easy to find via google if you want to see it).
  • vivek_
    PhD Student
    • Jul 2012
    • 164

    #2
    I assume the PL should be a log10 value, so a value of 0 indicates the most likely allele combination, which is BB In this case.

    SAMtools/BCFtools writes genotype likelihoods in the PL format which is a comma delimited list of phred-scaled data likelihoods of each possible genotype. For example, suppose REF=C and ALT=A,G, PL=7,0,37,13,40,49 means for the sample we are looking at, 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 and P(D|GG)=10^{-4.9}. This ordering has been changed since r921.
    Here's the info from the samtools mpileup page.

    Last edited by vivek_; 09-17-2012, 07:32 AM.

    Comment

    • swbarnes2
      Senior Member
      • May 2008
      • 910

      #3
      If your genotype ahs a quality of 3, you can't believe it, so yeah, treating that point as uncalled would make sense.

      The PLs are being given correctly. In the PLs a lower number is better quality. THE QUAL score in the QUAL column is backwards; that number represents the odds that the call is wrong, so higher is better there.

      Comment

      • jmartin
        Member
        • Dec 2009
        • 78

        #4
        Thank you for the answers, I feel like a dope forgetting that the PL values were log10 values

        I was wondering if there was a rule of thumb for what Genotype Quality (GQ) score people would consider 'good'?


        Thanks,
        John

        Comment

        Latest Articles

        Collapse

        ad_right_rmr

        Collapse

        News

        Collapse

        Topics Statistics Last Post
        Started by SEQadmin2, Yesterday, 10:09 AM
        0 responses
        10 views
        0 reactions
        Last Post SEQadmin2  
        Started by SEQadmin2, 06-04-2026, 08:59 AM
        0 responses
        19 views
        0 reactions
        Last Post SEQadmin2  
        Started by SEQadmin2, 06-02-2026, 12:03 PM
        0 responses
        26 views
        0 reactions
        Last Post SEQadmin2  
        Started by SEQadmin2, 06-02-2026, 11:40 AM
        0 responses
        21 views
        0 reactions
        Last Post SEQadmin2  
        Working...