SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
How to deal with no calls from GATK Unifiedgenotyper for indels audqf Bioinformatics 2 02-01-2012 02:53 PM
Using UnifiedGenotyper to call variants from haploid data? oiiio Bioinformatics 0 01-24-2012 08:29 PM
GATK update changes UnifiedGenotyper output doc.ramses Bioinformatics 2 10-18-2011 06:41 AM
GATK UnifiedGenotyper calling way too many SNPs in vcf swbarnes2 Bioinformatics 0 08-17-2011 01:33 PM
GATK UnifiedGenotyper SNP calling with Agilent 50Mb target enrichment nullalleles Genomic Resequencing 1 07-18-2011 06:28 AM

Reply
 
Thread Tools
Old 10-27-2011, 05:16 AM   #1
krobasky
Junior Member
 
Location: Boston, MA

Join Date: Dec 2010
Posts: 3
Default GATK UnifiedGenotyper not calling any variants

I'm trying to align RNASeq data to hg19 and call variants. Qualities look ok, the tophat alignment has been sorted (picard) and tagged with readgroups as per GATK discussions, and it all looks fine in IGV.

However not a single variant is being called. Here's the message that's most indicative of the problem (confidently called bases = 0):

Code:
INFO  23:00:50,322 UnifiedGenotyper - Visited bases                                3101976562
INFO  23:00:50,323 UnifiedGenotyper - Callable bases                               2864957043
INFO  23:00:50,323 UnifiedGenotyper - Confidently called bases                     0
INFO  23:00:50,323 UnifiedGenotyper - % callable bases of all loci                 92.359
INFO  23:00:50,323 UnifiedGenotyper - % confidently called bases of all loci       0.000
INFO  23:00:50,324 UnifiedGenotyper - % confidently called bases of callable loci  0.000
INFO  23:00:50,324 UnifiedGenotyper - Actual calls made                            0
INFO  23:00:50,325 TraversalEngine - Total runtime 14151.69 secs, 235.86 min, 3.93 hours
INFO  23:00:50,325 TraversalEngine - 66936259 reads were filtered out during traversal out of 77587056 total (86.27%)
INFO  23:00:50,326 TraversalEngine -   -> 66936259 reads (86.27% of total) failing MappingQualityUnavailableReadFilter
Would this kind of message seem to indicate a potential problem with the hg19 reference that I'm using?

I'm using:
Code:
ftp://ftp.broadinstitute.org/pub/seq/references/Homo_sapiens_assembly19.fasta
I call GATK as follows:
Code:
java -jar GenomeAnalysisTK.jar -I karyotypicRG.bam \
   -R Homo_sapiens_assembly19.fasta -T UnifiedGenotyper \
    -o snpCalls.vcf  \
    -stand_call_conf 50.0 \
    -stand_emit_conf 1 \
    -dcov 5000 >& gatkGenotyper.out
Adding a dbSNP ROD doesn't fix the problem, nor does EMIT_ALL_SITES (e.g., ignore quality scores and call everything).

Thanks in advance for any thoughts.
krobasky is offline   Reply With Quote
Old 10-27-2011, 08:41 AM   #2
raonyguimaraes
Member
 
Location: Belo Horizonte - Brazil

Join Date: Jun 2010
Posts: 38
Default

May be there is something wrong with the score of your FASTQ files. Are you using Illumina or Phred Scores ? You could try to use BWA for the alignment.
raonyguimaraes is offline   Reply With Quote
Old 10-28-2011, 02:07 PM   #3
krobasky
Junior Member
 
Location: Boston, MA

Join Date: Dec 2010
Posts: 3
Default

I'm using Illumina scores and tophat to find splicing junctions, and tophat only used bowtie, but I think this is a file format issue, somehow related to the reference.
krobasky is offline   Reply With Quote
Old 10-29-2011, 05:42 AM   #4
raonyguimaraes
Member
 
Location: Belo Horizonte - Brazil

Join Date: Jun 2010
Posts: 38
Default

First of all, extract metrics from your Alignment to check if they are ok, if they are may be it's because you are using a very high coverage 5000 means 500X of coverage ... Are you sure you have this ? Try to reduce to 300 ... and stand_emit_conf should not be 1 ... this is wrong ... use at least 10
raonyguimaraes is offline   Reply With Quote
Old 10-30-2011, 10:17 AM   #5
krobasky
Junior Member
 
Location: Boston, MA

Join Date: Dec 2010
Posts: 3
Default

I solved my problem by backing up to an earlier GATK version. Here's my hypothesis for what caused the problem:

After upgrading GATK, the UnifiedGenotyper began filtering most of my tophat-assembled reads. Spot-checking the alignment for quality scores (thanks for the tip, raonyguimaraes), I discovered that many of the higher-quality alignments received 255 for a mapping score (MAPQ), despite the following two assertions made in the SAM spec:

"No alignments should be assigned mapping quality 255."
"A value 255 indicates that the mapping quality is not available."

My guess is that an older version of GATK agrees with tophat's interpretation of 255, but more current versions follow the spec more closely.

I haven't completely chased down this hypothesis, but I did back up to an earlier version of GATK (v1.0.5974, Compiled 2011/06/10 13:26:59), and I no longer have problems with variant-calling my tophat alignment.

Thanks for the help!
krobasky is offline   Reply With Quote
Old 10-30-2011, 12:40 PM   #6
cjp
Member
 
Location: Cambridge, United Kingdom

Join Date: Jun 2011
Posts: 58
Default

The same thing happened for me with the latest GATK, so I had to use an old version to call SNPs in TopHat alignments.

For release 1.1 of GATK, there is something here about filtering out mappings with a read quality of 255.

http://www.broadinstitute.org/gsa/wi...TK_release_1.1

Chris
cjp is offline   Reply With Quote
Old 06-15-2012, 08:27 AM   #7
golharam
Member
 
Location: Philadelphia, PA

Join Date: Dec 2009
Posts: 55
Default any change to GATK for this?

I recently came across this problem with GATK 1.4-14-g2e47336. I assume this was changed a long time ago. Tophat still output 255 for mapping quality. Does anyone know if tophat or GATK has changed to accommodate this?
golharam is offline   Reply With Quote
Old 10-25-2012, 10:19 PM   #8
zhangfan
Junior Member
 
Location: china

Join Date: Apr 2010
Posts: 4
Default

Have you found the problem, I came across the same situation, everything looks fine, except for the "Actual calls made 0"
zhangfan 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 09:31 AM.


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