Bug Bulletin: The recent 3.2 release fixes many issues. If you run into a problem, please try the latest version before posting a bug report, as your problem may already have been solved.

Empty VCF file from GATK-Lite UnifiedGenotyper

rpandyarpandya Posts: 4Member

(There was another question about a similar symptom, but the answer does not appear to apply to what I'm seeing.)

I get an empty VCF file that just contains the header lines. The input VCF file is 1.8Gb, and as far as I can tell the content is OK - it has MAPQ scores, the flags seem reasonable, etc. I've attached a copy of the console output, and the beginning of the input file in SAM format. Let me know if you have any suggestions. Thanks,

Ravi

txt
txt
chr20-1204.console.txt
4K
txt
txt
out_realignment_chr20-1204.sam.txt
34K

Best Answer

Answers

  • Geraldine_VdAuweraGeraldine_VdAuwera Posts: 5,984Administrator, GATK Developer admin

    Hi Ravi,

    The first thing to try is to validate the input BAM file (using Picard's validateSamFile), just in case something's wrong with it. You can also look at the data in the IGV genome browser to visually check that the distribution of the reads looks OK. Then, try running your command again without -nt. Some filesystems don't do well with multithreading.

    Geraldine Van der Auwera, PhD

  • rpandyarpandya Posts: 4Member
    edited December 2012

    Thanks for the suggestions:

    • the reads look OK in IGV
    • removing -nt gives the same empty VCF file
    • Picard validateSamFile does complain about missing mate pairs after ~10M reads:

      C:\ravip\data>"\Program Files\java\jre7\bin\java.exe" -jar ..\bin\picard\ValidateSamFile.jar I=out_realignment_chr20-120 4.bam [Fri Dec 14 12:48:47 PST 2012] net.sf.picard.sam.ValidateSamFile INPUT=out_realignment_chr20-1204.bam MODE=VERBOSE MA X_OUTPUT=100 IGNORE_WARNINGS=false VALIDATE_INDEX=true IS_BISULFITE_SEQUENCED=false MAX_OPEN_TEMP_FILES=8000 VERBOSITY=I NFO QUIET=false VALIDATION_STRINGENCY=STRICT COMPRESSION_LEVEL=5 MAX_RECORDS_IN_RAM=500000 CREATE_INDEX=false CREATE_MD5 FILE=false [Fri Dec 14 12:48:47 PST 2012] Executing as ravip@ENDURE34 on Windows Server 2012 6.2 amd64; Java HotSpot(TM) 64-Bit Ser ver VM 1.7.0_09-b05; Picard version: 1.80(1295) INFO 2012-12-14 12:49:19 SamFileValidator Validated Read 10,000,000 records. Elapsed time: 00:00:31s. Time for last 10,000,000: 31s. Last read position: chr20:37,511,601 ERROR: Read name chr20_60726852_60727458?:?:??:?:?_MATERNAL_60764569_60765175_9ebbf, Mate not found for paired read ERROR: Read name chr20_11098955_11099486?:?:??:?:?_MATERNAL_11094073_11094604_2c2b9, Mate not found for paired read ERROR: Read name chr20_17146435_17147014?:?:??:?:?_MATERNAL_17143042_17143621_124911, Mate not found for paired read ERROR: Read name chr20_54530135_54530821?:?:??:?:?_PATERNAL_54562054_54562740_80fa24, Mate not found for paired read ERROR: Read name chr20_46667744_46668247?:?:??:?:?_MATERNAL_46662917_46663428_33771e, Mate not found for paired read ERROR: Read name chr20_1269711_1270225?:?:?_?:?:?_MATERNAL_1269623_1270117_2e7987, Mate not found for paired read .... Maximum output of [100] errors reached.

    Post edited by Geraldine_VdAuwera on
  • Geraldine_VdAuweraGeraldine_VdAuwera Posts: 5,984Administrator, GATK Developer admin

    Hi Ravi,

    Are you running GATK on Windows? Unfortunately we don't support that; we know that it's possible to get some tools to run, but there may be unexpected issues, and we have no way to help you troubleshoot them.

    Geraldine Van der Auwera, PhD

  • rpandyarpandya Posts: 4Member

    Hi Geraldine,

    Yes, I am running on Windows. However, I did try running the same commands on a Linux box and got the same results. Is the mate pair error a problem for UnifiedGenotyper? There's an earlier step of filtering reads that would be causing this that I can try removing.

  • Geraldine_VdAuweraGeraldine_VdAuwera Posts: 5,984Administrator, GATK Developer admin

    The mate pair error shouldn't be causing this; the UG applies a BadMate read filter that should take care of them. Are you running the latest version?

    Geraldine Van der Auwera, PhD

  • rpandyarpandya Posts: 4Member

    I just downloaded it: v2.2-10-gd1f076f

  • Geraldine_VdAuweraGeraldine_VdAuwera Posts: 5,984Administrator, GATK Developer admin

    We're actually on 2.2-16 right now, but that's close enough. We're in the process of releasing 2.3 so you might as well wait until it's out to update.

    Can you try calling variants on another file (maybe use a file from our resource bundle) to check if the UG works at all on your setup?

    Geraldine Van der Auwera, PhD

  • jiaqijiaqi Posts: 1

    Hi, I have met the same problem, I get an empty VCF file that just contains the header lines. Do you have any suggestion? Thansk a lot!

  • Geraldine_VdAuweraGeraldine_VdAuwera Posts: 5,984Administrator, GATK Developer admin

    Hi @jiaqi, can you please post your command line?

    Geraldine Van der Auwera, PhD

  • flyingflyersflyingflyers Posts: 13Member
    edited January 2013

    Me too!

    The Genome Analysis Toolkit (GATK) v2.0-38-g45f7b0d, Compiled 2012/08/09 22:00:16
    
    java -jar GenomeAnalysisTK.jar -T UnifiedGenotyper -R ../ref/GATK_PicardIndex/hg19c.fa -I bt2.fxmt.srt.dup.realign.bam -o bt2.fxmt.srt.dup.realign.dbsnp137.vcf -D ../ref/dbSNP/dbsnp_137.hg19.vcf -L ../ref/BEDfile/ExonOnly_50MbWE_ONTARGET.bed -glm SNP -mbq 30 -out_mode EMIT_VARIANTS_ONLY -gt_mode GENOTYPE_GIVEN_ALLELES
    

    output on screen

    ...
    INFO  19:15:51,345 ArgumentTypeDescriptor - Dynamically determined type of ../ref/dbSNP/dbsnp_137.hg19.vcf to be VCF
    INFO  19:15:51,356 GenomeAnalysisEngine - Strictness is SILENT
    INFO  19:15:51,429 SAMDataSource$SAMReaders - Initializing SAMRecords in serial
    INFO  19:15:51,453 SAMDataSource$SAMReaders - Done initializing BAM readers: total time 0.02
    INFO  19:15:51,472 RMDTrackBuilder - Loading Tribble index from disk for file ../ref/dbSNP/dbsnp_137.hg19.vcf
    INFO  19:15:54,061 TraversalEngine - [INITIALIZATION COMPLETE; TRAVERSAL STARTING]
    INFO  19:15:54,062 TraversalEngine -        Location processed.sites  runtime per.1M.sites completed total.runtime remaining
    INFO  19:16:23,937 TraversalEngine -   chr1:28081767        3.09e+06   30.1 s        9.7 s      1.9%        25.7 m    25.2 m
    ...
    INFO  19:42:55,666 TraversalEngine -  chrX:134855437        1.27e+08   27.0 m       12.7 s     99.1%        27.3 m    15.0 s
    INFO  19:43:08,405 TraversalEngine - Total runtime 1634.54 secs, 27.24 min, 0.45 hours
    INFO  19:43:08,405 TraversalEngine - 28244057 reads were filtered out during traversal out of 210793838 total (13.40%)
    INFO  19:43:08,406 TraversalEngine -   -> 1531326 reads (0.73% of total) failing BadMateFilter
    INFO  19:43:08,406 TraversalEngine -   -> 26263913 reads (12.46% of total) failing DuplicateReadFilter
    INFO  19:43:08,406 TraversalEngine -   -> 448818 reads (0.21% of total) failing UnmappedReadFilter
    INFO  19:43:09,723 GATKRunReport - Uploaded run statistics report to AWS S3
    

    The vcf file looks like (I deleted ## because it will turn into black )

    fileformat=VCFv4.1
    FILTER=<ID=LowQual,Description="Low quality">
    ...
    UnifiedGenotyper="analysis_type=UnifiedGenotyper input_file=[bt2.fxmt.srt.dup.realign.bam] read_buffer_size=null phone_home=STANDARD gatk_key=null read_filter=[] intervals=[../ref/BEDfile/029720_D_BED_20101013_merged_ExonOnly_50MbWE_ONTARGET.bed] excludeIntervals=null interval_set_rule=UNION interval_merging=ALL interval_padding=0 reference_sequence=../ref/GATK_PicardIndex/hg19c.fa nonDeterministicRandomSeed=false downsampling_type=BY_SAMPLE downsample_to_fraction=null downsample_to_coverage=250 baq=OFF baqGapOpenPenalty=40.0 performanceLog=null useOriginalQualities=false BQSR=null quantize_quals=0 disable_indel_quals=false emit_original_quals=false preserve_qscores_less_than=6 defaultBaseQualities=-1 validation_strictness=SILENT unsafe=null num_threads=24 num_cpu_threads=null num_io_threads=null num_bam_file_handles=null read_group_black_list=null pedigree=[] pedigreeString=[] pedigreeValidationType=STRICT allow_intervals_with_unindexed_bam=false generateShadowBCF=false useSlowGenotypes=false repairVCFHeader=null logging_level=INFO log_to_file=null help=false genotype_likelihoods_model=SNP p_nonref_model=EXACT heterozygosity=0.0010 pcr_error_rate=1.0E-4 genotyping_mode=GENOTYPE_GIVEN_ALLELES output_mode=EMIT_VARIANTS_ONLY standard_min_confidence_threshold_for_calling=30.0 standard_min_confidence_threshold_for_emitting=30.0 noSLOD=false annotateNDA=false alleles=(RodBinding name= source=UNBOUND) min_base_quality_score=30 max_deletion_fraction=0.05 max_alternate_alleles=3 cap_max_alternate_alleles_for_indels=false min_indel_count_for_genotyping=5 min_indel_fraction_per_sample=0.25 indel_heterozygosity=1.25E-4 indelGapContinuationPenalty=10 indelGapOpenPenalty=45 indelHaplotypeSize=80 noBandedIndel=false indelDebug=false ignoreSNPAlleles=false allReadsSP=false ignoreLaneInfo=false reference_sample_calls=(RodBinding name= source=UNBOUND) reference_sample_name=null sample_ploidy=2 min_quality_score=1 max_quality_score=40 site_quality_prior=20 min_power_threshold_for_calling=0.95 min_reference_depth=100 exclude_filtered_reference_sites=false dbsnp=(RodBinding name=dbsnp source=../ref/dbSNP/dbsnp_137.hg19.vcf) comp=[] out=org.broadinstitute.sting.gatk.io.stubs.VariantContextWriterStub no_cmdline_in_header=org.broadinstitute.sting.gatk.io.stubs.VariantContextWriterStub sites_only=org.broadinstitute.sting.gatk.io.stubs.VariantContextWriterStub debug_file=null metrics_file=null annotation=[] excludeAnnotation=[] filter_mismatching_base_and_quals=false"
    contig=<ID=chr1,length=249250621,assembly=hg19>
    ...
    contig=<ID=chrX,length=155270560,assembly=hg19>
    contig=<ID=chrY,length=59373566,assembly=hg19>
    reference=file:///scratch0/igm2/fwong/1KG_Ctrl_Sample_100bp/MAPQ_Substitute/../ref/GATK_PicardIndex/hg19c.fa
    CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  NA10847_100.1.2
    (END)
    

    It only happened after I updated to GATK 2.0. Any clue? Thanks.

    Post edited by Geraldine_VdAuwera on
  • Geraldine_VdAuweraGeraldine_VdAuwera Posts: 5,984Administrator, GATK Developer admin

    @flyingflyers, please try with the latest version (2.3) and let us know if the problem persists.

    Geraldine Van der Auwera, PhD

  • baescbaesc Posts: 21Member

    Hi, I've got the same problem using the GATK 2.3... can anyone tell me what I'm doing wrong?

    java \ -Xmx15g \ -jar /opt/local/bin/GATK-2.3/GenomeAnalysisTK.jar \ -T HaplotypeCaller \ -R /Users/zws/zws/projekte/dataUngesichert/BAM/Januar_Fries/BTA24_ALL/Chr24N.fa \ -I /Users/zws/zws/projekte/dataUngesichert/BAM/Januar_Fries/BTA24_ALL/55069_chr24.bam \ -o output.vcf

    Cheers, Chris

  • Geraldine_VdAuweraGeraldine_VdAuwera Posts: 5,984Administrator, GATK Developer admin

    Hi @baesc, can you please try again with the newest version (2.5.2) and see if the problem persists?

    Geraldine Van der Auwera, PhD

  • baescbaesc Posts: 21Member
    edited May 2013

    Hi @Geraldine_VdAuwera, I've tried again with the 2.5.2 Version and now get an error message saying that my reference file doesn't match my .bam file (ERROR MESSAGE: Badly formed genome loc: Contig Chr1 given as location, but this contig isn't present in the Fasta sequence dictionary).... Wierd, samtools accepted the reference, and version 2.3 accepted the reference and .dict files as well... any suggestions? Cheers, Chris

    Post edited by baesc on
  • baescbaesc Posts: 21Member
    edited May 2013

    .... aha, I need to use the "full" reference (UMD31.fasta) and not the chromosomal one (Chr24.fa)... is that right? Its running now (with UMD31.fasta), but it started walking at Chr1. I'd only like to call Chr24, but I'll see what happens....

    Post edited by baesc on
  • Geraldine_VdAuweraGeraldine_VdAuwera Posts: 5,984Administrator, GATK Developer admin

    Hi Chris,

    That's right, you need to use the full reference, as you probably have the other contigs in your bam dictionary, even if you subsetted the data.

    What you can do to just call one contig is use -L Chr24.

    Geraldine Van der Auwera, PhD

  • alyzaalyza baltimorePosts: 5Member

    Hi, I am having a similar problem. I am running GATK2.6 to call SNPs and indels but the vcf file is empty with some header lines. I was able to call SNPs using samtools. Here is my command line:

    java -jar /data/ngsc2data_a6/NGSC_scripts/Exome/GenomeAnalysisTK-2.6-5-gba531bd/GenomeAnalysisTK.jar -R /data/ngsc2data_a6/NGSC_reference/b37/human_g1k_v37.fasta -T UnifiedGenotyper -I T4_realigned_fixed.bam -o T4_GATK.vcf -stand_call_conf 50.0 -stand_emit_conf 10.0 -dcov 500 -glm BOTH --dbsnp /data/ngsc2data_a6/NGSC_reference/b37/dbsnp_137.b37.vcf

    Any ideas on what I am doing wrong?

    Thanx, Alyza

  • Geraldine_VdAuweraGeraldine_VdAuwera Posts: 5,984Administrator, GATK Developer admin

    Hi Alyza,

    Can you post the console output or logs from the run? It may be that you have a lot of reads getting filtered out, and that would be logged in the output.

    Geraldine Van der Auwera, PhD

  • alyzaalyza baltimorePosts: 5Member

    Hi, Thanx for your response. I am attaching the console output.

    Thanx, Alyza

    txt
    txt
    GATK_console_output.txt
    59K
  • Geraldine_VdAuweraGeraldine_VdAuwera Posts: 5,984Administrator, GATK Developer admin

    Okay, it looks like your data's not getting filtered out, that's good. Have you tried validating your bam with Picard's validateSamFile tool?

    The next step would be to look at the output of the run you did with samtools, pick a few high-confidence SNPs, and then try to call just those sites with GATK, and use the EMIT_ALL output mode.

    Geraldine Van der Auwera, PhD

  • alyzaalyza baltimorePosts: 5Member
    edited August 2013

    Hi, I ran Picard's validateSamFile and it complains about "Mate not found for paired read" for about 100 reads.

    I then reran GATK for only 3 of the high confidence SNPs found through samtools with the EMIT_ALL output mode. I am attaching the console output.

    #CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  T4
    1   565006  .   T   .   .   .   .   GT  ./.
    1   565901  .   G   .   .   .   .   GT  ./.
    1   565976  .   C   .   .   .   .   GT  ./. 
    

    Thanx, Alyza

    txt
    txt
    GATK_high_confidence.txt
    5K
    Post edited by Geraldine_VdAuwera on
  • Geraldine_VdAuweraGeraldine_VdAuwera Posts: 5,984Administrator, GATK Developer admin

    Hi Alyza,

    According to the log, all your reads for that region are filtered out due to MappingQualityUnavailableFilter. I'm not sure why this didn't come up in the log of your previous run, but you should check that your reads all have a value for MAPQ.

    Geraldine Van der Auwera, PhD

  • alyzaalyza baltimorePosts: 5Member

    Hi, I am not sure why that region was filtered since every position has a MAPQ value.

    chr1 565006 . T C 81 . DP=48;VDB=0.0777;AF1=0.500;AFE=0.500;DP4=0,0,23,25;MQ=40 PL 255,144,255 chr1 565901 . G A 17.1 . DP=69;VDB=0.0777;AF1=0.500;AFE=0.490;DP4=0,0,38,31;MQ=39 PL 255,208,255 chr1 565976 . C T 26 . DP=66;VDB=0.0615;AF1=0.500;AFE=0.499;DP4=0,0,21,45;MQ=40 PL 255,199,255

    Thanx, Alyza

  • Geraldine_VdAuweraGeraldine_VdAuwera Posts: 5,984Administrator, GATK Developer admin

    Hi Alyza,

    What you posted here are the variant records; what concerns me is that it looks like your reads don't have MAPQ. Can you please post a few entries from your SAM file?

    Geraldine Van der Auwera, PhD

  • alyzaalyza baltimorePosts: 5Member

    Hi.

    Below is a small subset of my sam file. The alignment was done using bowtie.

    HWI-ST1191:101:C2638ACXX:1:1101:1051:2050 163 5 35954223 255 50M = 35954381 208 AGAGGGGTGGCGTGTGCTGGGGTGGGGAGAACCTTCAAAGGACCAGGGAT =@@;BAD0<:F?C?EDHBB@GI68?@;FG8;36=?ACHCD?;@A;53(5 XA:i:0 MD:Z:50 NM:i:0 HWI-ST1191:101:C2638ACXX:1:1101:1051:2050 83 5 35954381 255 50M = 35954223 -208 TGCCCAGAGTGAGCCCCAGCAGAAACACAAAGACATCAATGAGGTACTGN B?>DB9?4FD;C6@IFGHCDD?C<3GC<>BE9FC<?@FDADDDBB=:1# XA:i:1 MD:Z:49C0 NM:i:1 HWI-ST1191:101:C2638ACXX:1:1101:1195:2092 99 11 66134954 255 50M = 66135034 130 TTTGGTCTCCAGCTCCTGAGCTTGGGCCTGGGATGATTTATTGGCCAGGT @@@BDBDFGHHHHJJGIIHGGIFHHIGIJJJJCHHIHIIGIJJIEGIIIE XA:i:0 MD:Z:50 NM:i:0 HWI-ST1191:101:C2638ACXX:1:1101:1195:2092 147 11 66135034 255 50M = 66134954 -130 GAAAGTGTTGGGCCTGCCTGGCTTGCCTGAGGTGTCTGGATCCCTCCCCA GHEGC@FC@@F@GB)BB;EIGGAHBHEFIGEHHD>IIHHAFHFDBFFC@@ XA:i:0 MD:Z:50 NM:i:0 HWI-ST1191:101:C2638ACXX:1:1101:1248:2155 99 8 86787860 255 50M = 86787958 148 CGAGGCTCTCCTTGCGGCCGTCCCGCACGTGCTGCTTTGCCACCTGGCAG CCCFFF

    Thanx, Alyza

Sign In or Register to comment.