SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
Interpreting Pindel output bwubb Bioinformatics 11 07-07-2014 05:47 AM
VarScan output to VCF OlgaMikh Bioinformatics 4 06-07-2012 11:46 PM
interpreting cuffdiff output with and without replicates pageskipro Bioinformatics 3 11-19-2011 12:25 AM
SNP vcf output: GT vs PL shuang Bioinformatics 0 08-09-2011 02:08 PM
Interpreting Genomica Module Networks output AbrahamB Bioinformatics 0 02-16-2011 11:11 AM

Reply
 
Thread Tools
Old 09-16-2012, 04:57 PM   #1
jmartin
Member
 
Location: St. Louis

Join Date: Dec 2009
Posts: 74
Default 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).
jmartin is offline   Reply With Quote
Old 09-17-2012, 07:22 AM   #2
vivek_
PhD Student
 
Location: Denmark

Join Date: Jul 2012
Posts: 164
Default

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.

Quote:
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.

http://samtools.sourceforge.net/mpileup.shtml

Last edited by vivek_; 09-17-2012 at 07:32 AM.
vivek_ is offline   Reply With Quote
Old 09-17-2012, 08:28 AM   #3
swbarnes2
Senior Member
 
Location: San Diego

Join Date: May 2008
Posts: 912
Default

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.
swbarnes2 is offline   Reply With Quote
Old 09-18-2012, 10:26 AM   #4
jmartin
Member
 
Location: St. Louis

Join Date: Dec 2009
Posts: 74
Default

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
jmartin is offline   Reply With Quote
Reply

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 04:52 PM.


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