SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
error during GATK indel realigner david.tamborero Bioinformatics 6 07-18-2012 06:30 AM
problem about GATK indel VQSR wanguan2000 Bioinformatics 2 11-07-2011 07:15 AM
GATK indel settings help please inou13 Bioinformatics 2 10-25-2011 09:09 AM
Good Indel Calling Options? Hkins552 Bioinformatics 1 08-24-2011 01:33 PM
SNP and indel calling methods Rachelly Bioinformatics 1 11-24-2010 04:45 AM

Reply
 
Thread Tools
Old 11-18-2011, 10:41 AM   #1
Alex Renwick
Member
 
Location: Houston, Texas

Join Date: Jul 2011
Posts: 44
Default gatk indel calling problem

I'm trying to get variant calling to work using a small set of simulated reads of the phi X genome. The problem I'm seeing is that indels are assigned a quality of zero and their vcf entries are incorrect.

I used wgsim from samtools to generate simulated paired-end fastq read files of varying depths of coverage, then used bwa to align the simulated reads to the phi X reference. I followed the GATK recommendations to realign around possible indels, then to recalibrate, then I ran the UnifiedGenotyper using the -glm INDEL option along with --output-mode EMIT_ALL_SITES.

The analysis info shows something like this:

Quote:
INFO 12:20:08,764 UnifiedGenotyper - Visited bases 5416
INFO 12:20:08,765 UnifiedGenotyper - Callable bases 7
INFO 12:20:08,765 UnifiedGenotyper - Confidently called bases 0
INFO 12:20:08,765 UnifiedGenotyper - % callable bases of all loci 0.129
INFO 12:20:08,765 UnifiedGenotyper - % confidently called bases of all loci 0.000
INFO 12:20:08,766 UnifiedGenotyper - % confidently called bases of callable loci 0.000
INFO 12:20:08,766 UnifiedGenotyper - Actual calls made 7
INFO 12:20:08,766 TraversalEngine - Total runtime 1.08 secs, 0.02 min, 0.00 hours

The corresponding vcf file shows (I added a space after the "GT:" to keep from happening):

Quote:
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT somesample
phix.fa 869 . T . 0 LowQual DP=183;MQ=48.18;MQ0=0;QD=0.00 GT: DP:PL ./.:183:0,0,0
phix.fa 875 . GCT . 0 LowQual DP=190;MQ=49.29;MQ0=0;QD=0.00 GT: DP:PL ./.:190:0,0,0
phix.fa 1408 . T . 0 LowQual DP=209;MQ=52.38;MQ0=0;QD=0.00 GT: DP:PL ./.:209:0,0,0
phix.fa 1409 . A . 0 LowQual DP=213;MQ=52.17;MQ0=0;QD=0.00 GT: DP:PL ./.:213:0,0,0
phix.fa 1949 . G . 0 LowQual DP=205;MQ=57.70;MQ0=0;QD=0.00 GT: DP:PL ./.:205:0,0,0
phix.fa 1950 . A . 0 LowQual DP=205;MQ=57.82;MQ0=0;QD=0.00 GT: DP:PL ./.:205:0,0,0
phix.fa 4225 . AC . 0 LowQual DP=179;MQ=57.51;MQ0=0;QD=0.00 GT: DP:PL ./.:179:0,0,0
The vcf entries correspond to actual indels in the simulated data. For example, at position 870 there is an insertion of CC, and 876 there is a deletion of CT. The problem is the the "ALT" column shows a "." where it should have the actual alternate allele (TCC and G in those two cases).

What can I do to 1) make the quality not zero; and 2) get the vcf file to correctly report the alternative allele?
Alex Renwick is offline   Reply With Quote
Old 11-18-2011, 11:16 AM   #2
NGSfan
Senior Member
 
Location: Austria

Join Date: Apr 2009
Posts: 181
Default

Quote:
Originally Posted by Alex Renwick View Post
I'm trying to get variant calling to work using a small set of simulated reads of the phi X genome. The problem I'm seeing is that indels are assigned a quality of zero and their vcf entries are incorrect.

I used wgsim from samtools to generate simulated paired-end fastq read files of varying depths of coverage, then used bwa to align the simulated reads to the phi X reference. I followed the GATK recommendations to realign around possible indels, then to recalibrate, then I ran the UnifiedGenotyper using the -glm INDEL option along with --output-mode EMIT_ALL_SITES.

The analysis info shows something like this:




The corresponding vcf file shows (I added a space after the "GT:" to keep from happening):



The vcf entries correspond to actual indels in the simulated data. For example, at position 870 there is an insertion of CC, and 876 there is a deletion of CT. The problem is the the "ALT" column shows a "." where it should have the actual alternate allele (TCC and G in those two cases).

What can I do to 1) make the quality not zero; and 2) get the vcf file to correctly report the alternative allele?
Hmmm.. you want to call all variants regardless of quality? try removing EMIT ALL SITES and lower the stand_emit_conf to 1
NGSfan is offline   Reply With Quote
Old 11-18-2011, 12:44 PM   #3
Alex Renwick
Member
 
Location: Houston, Texas

Join Date: Jul 2011
Posts: 44
Default

Quote:
Originally Posted by NGSfan View Post
Hmmm.. you want to call all variants regardless of quality? try removing EMIT ALL SITES and lower the stand_emit_conf to 1
I tried that--lowering the standard to zero, in fact--and without EMIT_ALL_SITES there is no output. It's like the indels are treated as nonvariant sites.
Alex Renwick is offline   Reply With Quote
Old 11-18-2011, 03:08 PM   #4
oiiio
Senior Member
 
Location: USA

Join Date: Jan 2011
Posts: 105
Default

Hi guys,
I have this exact problem. And I have tried -mmq 0 -mbq 0 -stand_emit_conf 0 -stand_call_conf 0 without seeing anything. I welcome any verification of my (non) findings with these parameters.

I also tried using real base qualities from real data without seeing any change... still actively working on this, and I would like to hear any solution or suggestion people have.

Last edited by oiiio; 11-18-2011 at 03:33 PM.
oiiio is offline   Reply With Quote
Old 12-02-2011, 11:56 AM   #5
Alex Renwick
Member
 
Location: Houston, Texas

Join Date: Jul 2011
Posts: 44
Default

I've done some experimenting and discovered that the problem comes from the amount of sequencing error in the input data.

The data I am using is simulated using wgsim (distributed with samtools). The default error rate is .02. With that rate, UnifiedGenotyper makes no indel calls, irrespective of setting call / emit confidences to zero and minimum base quality to zero. Snp calls work just fine, though.

UnifiedGenotyper is happy to make indel calls when the error rate is set to zero, and it works about half the time with an error rate of 0.01.

I haven't found any parameters that allow the user to specify sensitivity to sequencing errors. Anybody else have ideas?
Alex Renwick is offline   Reply With Quote
Old 12-02-2011, 12:17 PM   #6
oiiio
Senior Member
 
Location: USA

Join Date: Jan 2011
Posts: 105
Default

Nice find Alex,

When you say it works half the time, do you mean that it is calling about half the indels in the sample? Or that half the samples you run UnifiedGenotyper on are making calls at all?
oiiio is offline   Reply With Quote
Old 12-02-2011, 12:25 PM   #7
Alex Renwick
Member
 
Location: Houston, Texas

Join Date: Jul 2011
Posts: 44
Default

I've done around a dozen trials, and it seems to be all or nothing. Either all candidate sites are called, or none of them are.

edit: An error rate of 0.005 seems ok, too.

Last edited by Alex Renwick; 12-02-2011 at 12:28 PM.
Alex Renwick is offline   Reply With Quote
Reply

Tags
gatk, indel, unifiedgenotyper

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 11:24 PM.


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