Using Queue to run HaplotypeCaller and UnifiedGenotyper via shell script

rfarrerrfarrer Cambridge MAMember
edited November 2013 in Ask the GATK team

Hi,
I am trying to test the 2 SNP-callers of GATK (HC and UG) before choosing the right tool to use for the rest of my isolates. I have aligned my Illumina reads using BWA mem v.0.7.4-r385 (which is giving some unusual alignments with CIGAR codes beginning with deletions, which is causing various issues with other tools, but seems to be dealt with by -BadCigar in GATK, and is perhaps a separate issue). I am calling SNPs with GenomeAnalysisTK-2.7-4-g6f46d11, and I obtained Queue using git clone git://github.com/broadgsa/gatk-protected Sting, cd Sting, ant queue.

My shell script appears to work fine for the initial variant calling (for both HC and UG) before BQSR, but not after.
For UG, I get the following sort of output:

supercont_1.1 1 . A . . . DP=18;MQ=60.00;MQ0=0 GT ./.

supercont_1.1 2 . T . . . DP=28;MQ=60.00;MQ0=0 GT ./.

supercont_1.1 3 . T . . . DP=43;MQ=60.00;MQ0=0 GT ./.

supercont_1.1 4 . T . . . DP=53;MQ=60.00;MQ0=0 GT ./.

(I am using ref-base calling as well - but as you can see, no GT)
For HC, I do not even get an output, as it exits with this sort of error:

INFO 12:31:34,188 HelpFormatter - Date/Time: 2013/11/14 12:31:34
INFO 12:31:34,188 HelpFormatter - --------------------------------------------------------------------------------
INFO 12:31:34,188 HelpFormatter - --------------------------------------------------------------------------------
INFO 12:31:35,277 GenomeAnalysisEngine - Strictness is SILENT
INFO 12:31:35,449 GenomeAnalysisEngine - Downsampling Settings: Method: BY_SAMPLE, Target Coverage: 250
INFO 12:31:35,462 SAMDataSource$SAMReaders - Initializing SAMRecords in serial
INFO 12:31:35,504 SAMDataSource$SAMReaders - Done initializing BAM readers: total time 0.04
INFO 12:31:35,520 HCMappingQualityFilter - Filtering out reads with MAPQ < 20
INFO 12:31:35,538 IntervalUtils - Processing 35002 bp from intervals
INFO 12:31:35,768 GenomeAnalysisEngine - Preparing for traversal over 1 BAM files
INFO 12:31:35,890 GenomeAnalysisEngine - Done preparing for traversal
INFO 12:31:35,892 ProgressMeter - [INITIALIZATION COMPLETE; STARTING PROCESSING]
INFO 12:31:35,892 ProgressMeter - Location processed.active regions runtime per.1M.active regions completed total.runtime remaining
INFO 12:31:36,171 HaplotypeCaller - Using global mismapping rate of 45 => -4.5 in log10 likelihood units
INFO 12:32:05,921 ProgressMeter - supercont_1.1:491550 0.00e+00 30.0 s 49.6 w 4.3% 11.5 m 11.0 m
INFO 12:32:09,705 GATKRunReport - Uploaded run statistics report to AWS S3

ERROR ------------------------------------------------------------------------------------------
ERROR stack trace

java.lang.IllegalStateException: PairHMM Log Probability cannot be greater than 0: haplotype: [78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 67, 65, 71, 65, 71, 65, 84, 71, 65, 67, 71, 65, 65, 67, 67, 67, 67, 67, 71, 71, 65, 71, 65, 71, 65, 71, 67, 65, 71, 71, 71, 71, 65, 84, 65, 67, 67, 65, 67, 67, 67, 65, 71, 84, 67, 71, 84, 84, 65, 65, 71, 84, 65, 84, 67, 84, 65, 65, 84, 84, 67, 71, 65, 84, 71, 65, 84, 65, 84, 71, 67, 84, 67, 67, 65, 67, 65, 65, 84, 84, 67, 71, 65, 65, 84, 67, 71, 65, 67, 67, 84, 71, 71, 84, 84, 67, 71, 71, 71, 71, 84, 71, 84, 84, 71, 71, 84, 71, 71, 71, 65, 84, 65, 71, 65, 65, 65, 84, 84, 67, 84, 67, 71, 65, 65, 65, 84, 67, 67, 84, 84, 67, 84, 84, 65, 71, 71, 84, 67, 65, 65, 65, 65, 65, 65, 84, 65, 65, 65, 65, 65, 71, 84, 84, 71, 67, 71, 71, 84, 84, 71, 71, 71, 71, 84, 84, 65, 71, 65, 65, 84, 71, 71, 65, 71, 71, 67, 84, 67, 67, 84, 71, 67, 65, 65, 65, 67, 84, 71, 67, 84, 84, 71, 67, 65, 71, 65, 84, 71, 71, 71, 65, 71, 65, 65, 84, 71, 65, 67, 84, 84, 71, 71, 65, 84, 84, 65, 71, 84, 67, 67, 65, 71, 67, 84, 65, 84, 84, 84, 65, 84, 65, 84, 71, 71, 65, 71, 71, 71, 84, 65, 84, 67, 67, 71, 67, 84, 67, 71, 67, 71, 84], read: [71, 65, 84, 65, 84, 71, 67, 84, 67, 67, 65, 67, 65, 65, 84, 84, 67, 71, 65, 65, 84, 67, 71, 65, 67, 67, 84, 71, 71, 84, 84, 67, 71, 71, 71, 71, 84, 71, 84, 84, 84, 71, 71, 71, 84, 71, 65, 84, 65, 71, 65, 65, 65, 84, 84, 67, 84, 67, 71, 65, 65, 65, 84, 67, 67, 84, 84, 67, 84, 84, 65, 71, 71, 84, 67, 65, 65, 65, 65, 65, 65, 84, 65, 65, 65, 65, 65, 71, 84, 84, 71, 67, 71, 71, 84, 84, 71, 71, 71, 71], result: 0.002718
at org.broadinstitute.sting.utils.pairhmm.PairHMM.computeReadLikelihoodGivenHaplotypeLog10(PairHMM.java:131)
at org.broadinstitute.sting.gatk.walkers.haplotypecaller.LikelihoodCalculationEngine.computeReadLikelihoods(LikelihoodCalculationEngine.java:262)
at org.broadinstitute.sting.gatk.walkers.haplotypecaller.LikelihoodCalculationEngine.computeReadLikelihoods(LikelihoodCalculationEngine.java:205)
at org.broadinstitute.sting.gatk.walkers.haplotypecaller.HaplotypeCaller.map(HaplotypeCaller.java:770)
at org.broadinstitute.sting.gatk.walkers.haplotypecaller.HaplotypeCaller.map(HaplotypeCaller.java:140)
at org.broadinstitute.sting.gatk.traversals.TraverseActiveRegions$TraverseActiveRegionMap.apply(TraverseActiveRegions.java:708)
at org.broadinstitute.sting.gatk.traversals.TraverseActiveRegions$TraverseActiveRegionMap.apply(TraverseActiveRegions.java:704)
at org.broadinstitute.sting.utils.nanoScheduler.NanoScheduler.executeSingleThreaded(NanoScheduler.java:274)
at org.broadinstitute.sting.utils.nanoScheduler.NanoScheduler.execute(NanoScheduler.java:245)
at org.broadinstitute.sting.gatk.traversals.TraverseActiveRegions.traverse(TraverseActiveRegions.java:273)
at org.broadinstitute.sting.gatk.traversals.TraverseActiveRegions.traverse(TraverseActiveRegions.java:78)
at org.broadinstitute.sting.gatk.executive.LinearMicroScheduler.execute(LinearMicroScheduler.java:99)
at org.broadinstitute.sting.gatk.GenomeAnalysisEngine.execute(GenomeAnalysisEngine.java:313)
at org.broadinstitute.sting.gatk.CommandLineExecutable.execute(CommandLineExecutable.java:113)
at org.broadinstitute.sting.commandline.CommandLineProgram.start(CommandLineProgram.java:245)
at org.broadinstitute.sting.commandline.CommandLineProgram.start(CommandLineProgram.java:152)
at org.broadinstitute.sting.gatk.CommandLineGATK.main(CommandLineGATK.java:91)

ERROR ------------------------------------------------------------------------------------------
ERROR A GATK RUNTIME ERROR has occurred (version 2.7-4-g6f46d11):
ERROR
ERROR This might be a bug. Please check the documentation guide to see if this is a known problem.
ERROR If not, please post the error message, with stack trace, to the GATK forum.
ERROR Visit our website and forum for extensive documentation and answers to
ERROR commonly asked questions http://www.broadinstitute.org/gatk
ERROR
ERROR MESSAGE: PairHMM Log Probability cannot be greater than 0: haplotype: [78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 67, 65, 71, 65, 71, 65, 84, 71, 65, 67, 71, 65, 65, 67, 67, 67, 67, 67, 71, 71, 65, 71, 65, 71, 65, 71, 67, 65, 71, 71, 71, 71, 65, 84, 65, 67, 67, 65, 67, 67, 67, 65, 71, 84, 67, 71, 84, 84, 65, 65, 71, 84, 65, 84, 67, 84, 65, 65, 84, 84, 67, 71, 65, 84, 71, 65, 84, 65, 84, 71, 67, 84, 67, 67, 65, 67, 65, 65, 84, 84, 67, 71, 65, 65, 84, 67, 71, 65, 67, 67, 84, 71, 71, 84, 84, 67, 71, 71, 71, 71, 84, 71, 84, 84, 71, 71, 84, 71, 71, 71, 65, 84, 65, 71, 65, 65, 65, 84, 84, 67, 84, 67, 71, 65, 65, 65, 84, 67, 67, 84, 84, 67, 84, 84, 65, 71, 71, 84, 67, 65, 65, 65, 65, 65, 65, 84, 65, 65, 65, 65, 65, 71, 84, 84, 71, 67, 71, 71, 84, 84, 71, 71, 71, 71, 84, 84, 65, 71, 65, 65, 84, 71, 71, 65, 71, 71, 67, 84, 67, 67, 84, 71, 67, 65, 65, 65, 67, 84, 71, 67, 84, 84, 71, 67, 65, 71, 65, 84, 71, 71, 71, 65, 71, 65, 65, 84, 71, 65, 67, 84, 84, 71, 71, 65, 84, 84, 65, 71, 84, 67, 67, 65, 71, 67, 84, 65, 84, 84, 84, 65, 84, 65, 84, 71, 71, 65, 71, 71, 71, 84, 65, 84, 67, 67, 71, 67, 84, 67, 71, 67, 71, 84], read: [71, 65, 84, 65, 84, 71, 67, 84, 67, 67, 65, 67, 65, 65, 84, 84, 67, 71, 65, 65, 84, 67, 71, 65, 67, 67, 84, 71, 71, 84, 84, 67, 71, 71, 71, 71, 84, 71, 84, 84, 84, 71, 71, 71, 84, 71, 65, 84, 65, 71, 65, 65, 65, 84, 84, 67, 84, 67, 71, 65, 65, 65, 84, 67, 67, 84, 84, 67, 84, 84, 65, 71, 71, 84, 67, 65, 65, 65, 65, 65, 65, 84, 65, 65, 65, 65, 65, 71, 84, 84, 71, 67, 71, 71, 84, 84, 71, 71, 71, 71], result: 0.002718
ERROR ------------------------------------------------------------------------------------------

The relevant part of the shell script I have been using to generate these:

Qscript:

92 writeQscript() {

93 # variables: 1) program 2) out-VCF 3) out-QScript

94 QSCRIPT="package org.broadinstitute.sting.queue.qscripts.examples\n

95 import org.broadinstitute.sting.queue.QScript\n

96 import org.broadinstitute.sting.queue.extensions.gatk._\n

97 class Genotyper extends QScript {\n

98 qscript =>\n

99 // Reference & BAM\n

100 @Input(doc=\"Reference file\", shortName=\"R\")\n

101 var referenceFile: File = _\n

102 @Input(doc=\"Bam file\", shortName=\"I\")

103 var bamFile: File = _\n\n

104 @Argument(doc=\"Optional filter names.\", shortName=\"filter\", required=false)\n

105 var filterNames: List[String] ="

106 if [ $1 == "UnifiedGenotyper" ] ; then

107 QSCRIPT="$QSCRIPT List(\"output_mode=EMIT_ALL_SITES\", \"allSitePLs\")\n"

108 else

109 QSCRIPT="$QSCRIPT List(\"ERC\")\n"

110 fi

111 QSCRIPT="$QSCRIPT trait UnifiedGenotyperArguments extends CommandLineGATK {\n

112 this.reference_sequence = qscript.referenceFile\n

113 this.memoryLimit = 2 // For H.C. this was set to 12, for U.G., 2\n

114 }\n

115 def script() {\n

116 val genotyper = new $1 with UnifiedGenotyperArguments\n

117 genotyper.scatterCount = 500\n

118 genotyper.input_file :+= qscript.bamFile\n

119 genotyper.read_filter :+= \"BadCigar\"\n

120 genotyper.out = \"$2\"\n"

121 if [ $1 == "UnifiedGenotyper" ] ; then

122 QSCRIPT="$QSCRIPT genotyper.out_mode = org.broadinstitute.sting.gatk.walkers.genotyper.UnifiedGenotyperEngine.OUTPUT_MODE.EMIT_ALL_SITES\n"

123 else

124 QSCRIPT="$QSCRIPT genotyper.ERC = org.broadinstitute.sting.gatk.walkers.haplotypecaller.HaplotypeCaller.ReferenceConfidenceMode.GVCF\n"

125 fi

126 QSCRIPT="$QSCRIPT add(genotyper)\n}\n}"

127 echo -e $QSCRIPT > $3

128 }

129

130 if [ $3 == "all" ] || [ $3 == "bqsr" ] ; then

131 echo "GATK initial variant calling for BQSR inc. ref calling with $4"

132 writeQscript $4 $RAWVCF $QSCRIPTNAME

133 java -Djava.io.tmpdir=$TMPDIR -jar Queue.jar -S $QSCRIPTNAME -R $1 -I $BAMFOUR -bsub -jobQueue week -run

139 java -Xmx8g -jar GenomeAnalysisTK.jar -T BaseRecalibrator -R $1 -I $BAMFOUR -o $RECALIBRATION -knownSites $RAWVCF -nct 8

140 java -Xmx8g -jar GenomeAnalysisTK.jar -T PrintReads -R $1 -I $BAMFOUR -BQSR $RECALIBRATION -o $BAMFIVE -nct 4

141 fi

142

143 # Initial (currently final) variant discovery and genotyping

144 if [ $3 == "all" ] || [ $3 == "variant" ] ; then

145 echo "GATK Initial variant discovery and genotyping"

146 writeQscript $4 $OUTVCF $QSCRIPTNAME

147 java -Djava.io.tmpdir=$TMPDIR -jar Queue.jar -S $QSCRIPTNAME -R $1 -I $BAMFIVE -bsub -jobQueue week -run

148 fi

Any help with any of the issues I have been having would be great.
Thanks,
Rhys

Answers

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    Hi Rhys,

    For troubleshooting purposes, I would recommend you deconvolute your script into standalone command lines that you can test sequentially. This should help you pinpoint what step is problematic. Once you can show me a command line and corresponding results, and specifically indicate what is wrong with the results, I can help you fix it.

  • rfarrerrfarrer Cambridge MAMember

    Hi, thanks for the quick response.
    The specific command that is giving the problem for both HaplotypeCaller and UnifiedGenotyper is:

    java -Djava.io.tmpdir=$TMPDIR -jar Queue.jar -S $QSCRIPTNAME -R $1 -I $BAMFIVE -bsub -jobQueue week -run

    where:
    $1 is the reference fasta
    $TMPDIR is the bam file . '-tmp'
    $QSCRIPTNAME (for HaplotypeCaller here) made by the shell script is:

    package org.broadinstitute.sting.queue.qscripts.examples

    import org.broadinstitute.sting.queue.QScript

    import org.broadinstitute.sting.queue.extensions.gatk._

    class Genotyper extends QScript {

    qscript =>

    // Reference & BAM

    @Input(doc="Reference file", shortName="R")

    var referenceFile: File = _

    @Input(doc="Bam file", shortName="I") var bamFile: File = _

    @Argument(doc="Optional filter names.", shortName="filter", required=false)

    var filterNames: List[String] = List("ERC")

    trait UnifiedGenotyperArguments extends CommandLineGATK {

    this.reference_sequence = qscript.referenceFile

    this.memoryLimit = 2 // For H.C. this was set to 12, for U.G., 2

    }

    def script() {

    val genotyper = new HaplotypeCaller with UnifiedGenotyperArguments

    genotyper.scatterCount = 500

    genotyper.input_file :+= qscript.bamFile

    genotyper.read_filter :+= "BadCigar"

    genotyper.out = "Sample_HL_B8.sorted.sam_1.fastq-with-pairs-aln.sorted.bam.Post-BQSR-using-HaplotypeCaller.vcf"

    genotyper.ERC = org.broadinstitute.sting.gatk.walkers.haplotypecaller.HaplotypeCaller.ReferenceConfidenceMode.GVCF

    add(genotyper)

    } }

    For this command, it appears to split the jobs OK. Most run, but some do not. I have a very long output, but I think this is where the output starts to show errors:

    INFO 12:32:16,548 QGraph - 193 Pend, 309 Run, 0 Fail, 1 Done
    ERROR 12:32:17,767 FunctionEdge - Error: 'java' '-Xmx2048m' '-XX:+UseParallelOldGC' '-XX:ParallelGCThreads=4' '-XX:GCTimeLimit=50' '-XX:GCHeapFreeLimit=10' '-Djava.io.tmpdir=/seq/gscidC/organisms/Crypto_gattii/Imperial_Birmingham_sequences/Sample_HL_B8/FDR_60K_SNP/full/Sample_HL_B8.sorted.sam_1.fastq-with-pairs-aln.sorted.bam-new-headers-DupRem-reordered-realigned.bam-tmp' '-cp' '/seq/gscidC/organisms/Crypto_gattii/Imperial_Birmingham_sequences/Sample_HL_B8/FDR_60K_SNP/full/Queue.jar' 'org.broadinstitute.sting.gatk.CommandLineGATK' '-T' 'HaplotypeCaller' '-I' '/seq/gscidC/organisms/Crypto_gattii/Imperial_Birmingham_sequences/Sample_HL_B8/FDR_60K_SNP/full/Sample_HL_B8.sorted.sam_1.fastq-with-pairs-aln.sorted.bam-new-headers-DupRem-reordered-realigned-recal-HaplotypeCaller.bam' '-rf' 'BadCigar' '-L' '/seq/gscidC/organisms/Crypto_gattii/Imperial_Birmingham_sequences/Sample_HL_B8/FDR_60K_SNP/full/.queue/scatterGather/HaplotypeCaller-1-sg/temp_015_of_500/scatter.intervals' '-R' '/home/unix/rfarrer/Crypto_gattii/analysis-RF/R265_FDR_REFS/Cgattii_R265_genome.fa-60000-SNP-in-Genome.fasta' '-o' '/seq/gscidC/organisms/Crypto_gattii/Imperial_Birmingham_sequences/Sample_HL_B8/FDR_60K_SNP/full/.queue/scatterGather/HaplotypeCaller-1-sg/temp_015_of_500/Sample_HL_B8.sorted.sam_1.fastq-with-pairs-aln.sorted.bam.Post-BQSR-using-HaplotypeCaller.vcf' '-ERC' 'GVCF'
    ERROR 12:32:17,846 FunctionEdge - Contents of /seq/gscidC/organisms/Crypto_gattii/Imperial_Birmingham_sequences/Sample_HL_B8/FDR_60K_SNP/full/.queue/scatterGather/HaplotypeCaller-1-sg/temp_015_of_500/Sample_HL_B8.sorted.sam_1.fastq-with-pairs-aln.sorted.bam.Post-BQSR-using-HaplotypeCaller.vcf.out:
    Sender: LSF System <[email protected]>
    Subject: Job 4849861: < 'java' '-Xmx2048m' '-XX:+UseParallelOldGC' '-XX:ParallelGCThreads=4' '-XX:GCTimeLimit=50' '-XX:GCHeapFreeLimit=10' '-Djava.io.tmpdir=/seq/gscidC/organisms/Crypto_gattii/Imperial_Birmingham_sequences/Sample_HL_B8/FDR_60K_SNP/full/Sample_HL_B8.sorted.sam_1.fastq-with-pairs-aln.sorted.bam-new-headers-DupRem-reordered-realigned.bam-tmp' '-cp' '/seq/gscidC/organisms/Crypto_gattii/Imperial_Birmingham_sequences/Sample_HL_B8/FDR_60K_SNP/full/Queue.jar' 'org.broadinstitute.sting.gatk.CommandLineGATK' '-T' 'HaplotypeCaller' '-I' '/seq/gscidC/organisms/Crypto_gattii/Imperial_Birmingham_sequences/Sample_HL_B8/FDR_60K_SNP/full/Sample_HL_B8.sorted.sam_1.fastq-with-pairs-aln.sorted.bam-new-headers-DupRem-reordered-realigned-recal-HaplotypeCaller.bam' '-rf' 'BadCigar' '-L' '/seq/gscidC/organisms/Crypto_gattii/Imperial_Birmingham_sequences/Sample_HL_B8/FDR_60K_SNP/full/.queue/scatterGather/HaplotypeCaller-1-sg/temp_015_of_500/scatter.intervals' '-R' '/home/unix/rfarrer/Crypto_gattii/analysis-RF/R265_FDR_REFS/Cgattii_R265_genome.fa-60000-SNP-in-Genome.fasta' '-o' '/seq/gscidC/organisms/Crypto_gattii/Imperial_Birmingham_sequences/Sample_HL_B8/FDR_60K_SNP/full/.queue/scatterGather/HaplotypeCaller-1-sg/temp_015_of_500/Sample_HL_B8.sorted.sam_1.fastq-with-pairs-aln.sorted.bam.Post-BQSR-using-HaplotypeCaller.vcf' '-ERC' 'GVCF' > Exited

    Job < 'java' '-Xmx2048m' '-XX:+UseParallelOldGC' '-XX:ParallelGCThreads=4' '-XX:GCTimeLimit=50' '-XX:GCHeapFreeLimit=10' '-Djava.io.tmpdir=/seq/gscidC/organisms/Crypto_gattii/Imperial_Birmingham_sequences/Sample_HL_B8/FDR_60K_SNP/full/Sample_HL_B8.sorted.sam_1.fastq-with-pairs-aln.sorted.bam-new-headers-DupRem-reordered-realigned.bam-tmp' '-cp' '/seq/gscidC/organisms/Crypto_gattii/Imperial_Birmingham_sequences/Sample_HL_B8/FDR_60K_SNP/full/Queue.jar' 'org.broadinstitute.sting.gatk.CommandLineGATK' '-T' 'HaplotypeCaller' '-I' '/seq/gscidC/organisms/Crypto_gattii/Imperial_Birmingham_sequences/Sample_HL_B8/FDR_60K_SNP/full/Sample_HL_B8.sorted.sam_1.fastq-with-pairs-aln.sorted.bam-new-headers-DupRem-reordered-realigned-recal-HaplotypeCaller.bam' '-rf' 'BadCigar' '-L' '/seq/gscidC/organisms/Crypto_gattii/Imperial_Birmingham_sequences/Sample_HL_B8/FDR_60K_SNP/full/.queue/scatterGather/HaplotypeCaller-1-sg/temp_015_of_500/scatter.intervals' '-R' '/home/unix/rfarrer/Crypto_gattii/analysis-RF/R265_FDR_REFS/Cgattii_R265_genome.fa-60000-SNP-in-Genome.fasta' '-o' '/seq/gscidC/organisms/Crypto_gattii/Imperial_Birmingham_sequences/Sample_HL_B8/FDR_60K_SNP/full/.queue/scatterGather/HaplotypeCaller-1-sg/temp_015_of_500/Sample_HL_B8.sorted.sam_1.fastq-with-pairs-aln.sorted.bam.Post-BQSR-using-HaplotypeCaller.vcf' '-ERC' 'GVCF' > was submitted from host by user in cluster .
    Job was executed on host(s) , in queue , as user in cluster .
    </home/unix/rfarrer> was used as the home directory.
    </seq/gscidC/organisms/Crypto_gattii/Imperial_Birmingham_sequences/Sample_HL_B8/FDR_60K_SNP/full/.queue/scatterGather/HaplotypeCaller-1-sg/temp_015_of_500> was used as the working directory.
    Started at Thu Nov 14 12:31:23 2013
    Results reported at Thu Nov 14 12:32:10 2013

    Your job looked like:


    LSBATCH: User input

    sh /seq/gscidC/organisms/Crypto_gattii/Imperial_Birmingham_sequences/Sample_HL_B8/FDR_60K_SNP/full/Sample_HL_B8.sorted.sam_1.fastq-with-pairs-aln.sorted.bam-new-headers-DupRem-reordered-realigned.bam-tmp/.exec9117791809202266754

    Exited with exit code 1.

    Resource usage summary:

    CPU time   :     53.77 sec.
    Max Memory :       771 MB
    Max Swap   :      3313 MB
    
    Max Processes  :         4
    Max Threads    :        19
    

    The output (if any) follows:

    INFO 12:31:34,175 HelpFormatter - --------------------------------------------------------------------------------
    INFO 12:31:34,179 HelpFormatter - The Genome Analysis Toolkit (GATK) v2.7-4-g6f46d11, Compiled 2013/11/01 14:44:45
    INFO 12:31:34,179 HelpFormatter - Copyright (c) 2010 The Broad Institute
    INFO 12:31:34,180 HelpFormatter - For support and documentation go to http://www.broadinstitute.org/gatk
    INFO 12:31:34,187 HelpFormatter - Program Args: -T HaplotypeCaller -I /seq/gscidC/organisms/Crypto_gattii/Imperial_Birmingham_sequences/Sample_HL_B8/FDR_60K_SNP/full/Sample_HL_B8.sorted.sam_1.fastq-with-pairs-aln.sorted.bam-new-headers-DupRem-reordered-realigned-recal-HaplotypeCaller.bam -rf BadCigar -L /seq/gscidC/organisms/Crypto_gattii/Imperial_Birmingham_sequences/Sample_HL_B8/FDR_60K_SNP/full/.queue/scatterGather/HaplotypeCaller-1-sg/temp_015_of_500/scatter.intervals -R /home/unix/rfarrer/Crypto_gattii/analysis-RF/R265_FDR_REFS/Cgattii_R265_genome.fa-60000-SNP-in-Genome.fasta -o /seq/gscidC/organisms/Crypto_gattii/Imperial_Birmingham_sequences/Sample_HL_B8/FDR_60K_SNP/full/.queue/scatterGather/HaplotypeCaller-1-sg/temp_015_of_500/Sample_HL_B8.sorted.sam_1.fastq-with-pairs-aln.sorted.bam.Post-BQSR-using-HaplotypeCaller.vcf -ERC GVCF
    INFO 12:31:34,188 HelpFormatter - Date/Time: 2013/11/14 12:31:34
    INFO 12:31:34,188 HelpFormatter - --------------------------------------------------------------------------------
    INFO 12:31:34,188 HelpFormatter - --------------------------------------------------------------------------------
    INFO 12:31:35,277 GenomeAnalysisEngine - Strictness is SILENT
    INFO 12:31:35,449 GenomeAnalysisEngine - Downsampling Settings: Method: BY_SAMPLE, Target Coverage: 250
    INFO 12:31:35,462 SAMDataSource$SAMReaders - Initializing SAMRecords in serial
    INFO 12:31:35,504 SAMDataSource$SAMReaders - Done initializing BAM readers: total time 0.04
    INFO 12:31:35,520 HCMappingQualityFilter - Filtering out reads with MAPQ < 20
    INFO 12:31:35,538 IntervalUtils - Processing 35002 bp from intervals
    INFO 12:31:35,768 GenomeAnalysisEngine - Preparing for traversal over 1 BAM files
    INFO 12:31:35,890 GenomeAnalysisEngine - Done preparing for traversal
    INFO 12:31:35,892 ProgressMeter - [INITIALIZATION COMPLETE; STARTING PROCESSING]
    INFO 12:31:35,892 ProgressMeter - Location processed.active regions runtime per.1M.active regions completed total.runtime remaining
    INFO 12:31:36,171 HaplotypeCaller - Using global mismapping rate of 45 => -4.5 in log10 likelihood units
    INFO 12:32:05,921 ProgressMeter - supercont_1.1:491550 0.00e+00 30.0 s 49.6 w 4.3% 11.5 m 11.0 m
    INFO 12:32:09,705 GATKRunReport - Uploaded run statistics report to AWS S3

    ERROR ------------------------------------------------------------------------------------------
    ERROR stack trace

    java.lang.IllegalStateException: PairHMM Log Probability cannot be greater than 0: haplotype: [78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78,
    78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 67, 65, 71, 65, 71, 65, 84, 71, 65, 67, 71, 65, 65, 67, 67, 67, 67, 67, 71, 71, 65, 71, 65, 71, 65, 71, 67, 65, 71, 71, 71, 71, 65, 84, 65, 67, 67, 65, 67, 67, 67, 65, 71, 84, 67, 71, 84, 84, 65, 65, 71, 84, 65, 84, 67, 84, 65, 65, 84, 84, 67, 71, 65, 84, 71, 65, 84, 65, 84, 71, 67, 84, 67, 67, 65, 67, 65, 65, 84, 84, 67, 71, 65, 65, 84, 67, 71, 65, 67, 67, 84, 71, 71, 84, 84, 67, 71, 71, 71, 71, 84, 71, 84, 84, 71, 71, 84, 71, 71, 71, 65, 84, 65, 71, 65, 65, 65, 84, 84, 67, 84, 67, 71, 65, 65, 65, 84, 67, 67, 84, 84, 67, 84, 84, 65, 71, 71, 84, 67, 65, 65, 65, 65, 65, 65, 84, 65, 65, 65, 65, 65, 71, 84, 84, 71, 67, 71, 71, 84, 84, 71, 71, 71, 71, 84, 84, 65, 71, 65, 65, 84, 71, 71, 65, 71, 71, 67, 84, 67, 67, 84, 71, 67, 65, 65, 65, 67, 84, 71, 67, 84, 84, 71, 67, 65, 71, 65, 84, 71, 71, 71, 65, 71, 65, 65, 84, 71, 65, 67, 84, 84, 71, 71, 65, 84, 84, 65, 71, 84, 67, 67, 65, 71, 67, 84, 65, 84, 84, 84, 65, 84, 65, 84, 71, 71, 65, 71, 71, 71, 84, 65, 84, 67, 67, 71, 67, 84, 67, 71, 67, 71, 84], read: [71, 65, 84, 65, 84, 71, 67, 84, 67, 67, 65, 67, 65, 65, 84, 84, 67, 71, 65, 65, 84, 67, 71, 65, 67, 67, 84, 71, 71, 84, 84, 67, 71, 71, 71, 71, 84, 71, 84, 84, 84, 71, 71, 71, 84, 71, 65, 84, 65, 71, 65, 65, 65, 84, 84, 67, 84, 67, 71, 65, 65, 65, 84, 67, 67, 84, 84, 67, 84, 84, 65, 71, 71, 84, 67, 65, 65, 65, 65, 65, 65, 84, 65, 65, 65, 65, 65, 71, 84, 84, 71, 67, 71, 71, 84, 84, 71, 71, 71, 71], result: 0.002718
    at org.broadinstitute.sting.utils.pairhmm.PairHMM.computeReadLikelihoodGivenHaplotypeLog10(PairHMM.java:131)
    at org.broadinstitute.sting.gatk.walkers.haplotypecaller.LikelihoodCalculationEngine.computeReadLikelihoods(LikelihoodCalculationEngine.java:262)
    at org.broadinstitute.sting.gatk.walkers.haplotypecaller.LikelihoodCalculationEngine.computeReadLikelihoods(LikelihoodCalculationEngine.java:205)
    at org.broadinstitute.sting.gatk.walkers.haplotypecaller.HaplotypeCaller.map(HaplotypeCaller.java:770)
    at org.broadinstitute.sting.gatk.walkers.haplotypecaller.HaplotypeCaller.map(HaplotypeCaller.java:140)
    at org.broadinstitute.sting.gatk.traversals.TraverseActiveRegions$TraverseActiveRegionMap.apply(TraverseActiveRegions.java:708)
    at org.broadinstitute.sting.gatk.traversals.TraverseActiveRegions$TraverseActiveRegionMap.apply(TraverseActiveRegions.java:704)
    at org.broadinstitute.sting.utils.nanoScheduler.NanoScheduler.executeSingleThreaded(NanoScheduler.java:274)
    at org.broadinstitute.sting.utils.nanoScheduler.NanoScheduler.execute(NanoScheduler.java:245)
    at org.broadinstitute.sting.gatk.traversals.TraverseActiveRegions.traverse(TraverseActiveRegions.java:273)
    at org.broadinstitute.sting.gatk.traversals.TraverseActiveRegions.traverse(TraverseActiveRegions.java:78)
    at org.broadinstitute.sting.gatk.executive.LinearMicroScheduler.execute(LinearMicroScheduler.java:99)
    at org.broadinstitute.sting.gatk.GenomeAnalysisEngine.execute(GenomeAnalysisEngine.java:313)
    at org.broadinstitute.sting.gatk.CommandLineExecutable.execute(CommandLineExecutable.java:113)
    at org.broadinstitute.sting.commandline.CommandLineProgram.start(CommandLineProgram.java:245)
    at org.broadinstitute.sting.commandline.CommandLineProgram.start(CommandLineProgram.java:152)
    at org.broadinstitute.sting.gatk.CommandLineGATK.main(CommandLineGATK.java:91)

    ERROR ------------------------------------------------------------------------------------------
    ERROR A GATK RUNTIME ERROR has occurred (version 2.7-4-g6f46d11):
    ERROR
    ERROR This might be a bug. Please check the documentation guide to see if this is a known problem.
    ERROR If not, please post the error message, with stack trace, to the GATK forum.
    ERROR Visit our website and forum for extensive documentation and answers to
    ERROR commonly asked questions http://www.broadinstitute.org/gatk
    ERROR
    ERROR MESSAGE: PairHMM Log Probability cannot be greater than 0: haplotype: [78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 78, 67, 65, 71, 65, 71, 65, 84, 71, 65, 67, 71, 65, 65, 67, 67, 67, 67, 67, 71, 71, 65, 71, 65, 71, 65, 71, 67, 65, 71, 71, 71, 71, 65, 84, 65, 67, 67, 65, 67, 67, 67, 65, 71, 84, 67, 71, 84, 84, 65, 65, 71, 84, 65, 84, 67, 84, 65, 65, 84, 84, 67, 71, 65, 84, 71, 65, 84, 65, 84, 71, 67, 84, 67, 67, 65, 67, 65, 65, 84, 84, 67, 71, 65, 65, 84, 67, 71, 65, 67, 67, 84, 71, 71, 84, 84, 67, 71, 71, 71, 71, 84, 71, 84, 84, 71, 71, 84, 71, 71, 71, 65, 84, 65, 71, 65, 65, 65, 84, 84, 67, 84, 67, 71, 65, 65, 65, 84, 67, 67, 84, 84, 67, 84, 84, 65, 71, 71, 84, 67, 65, 65, 65, 65, 65, 65, 84, 65, 65, 65, 65, 65, 71, 84, 84, 71, 67, 71, 71, 84, 84, 71, 71, 71, 71, 84, 84, 65, 71, 65, 65, 84, 71, 71, 65, 71, 71, 67, 84, 67, 67, 84, 71, 67, 65, 65, 65, 67, 84, 71, 67, 84, 84, 71, 67, 65, 71, 65, 84, 71, 71, 71, 65, 71, 65, 65, 84, 71, 65, 67, 84, 84, 71, 71, 65, 84, 84, 65, 71, 84, 67, 67, 65, 71, 67, 84, 65, 84, 84, 84, 65, 84, 65, 84, 71, 71, 65, 71, 71, 71, 84, 65, 84, 67, 67, 71, 67, 84, 67, 71, 67, 71, 84], read: [71, 65, 84, 65, 84, 71, 67, 84, 67, 67, 65, 67, 65, 65, 84, 84, 67, 71, 65, 65, 84, 67, 71, 65, 67, 67, 84, 71, 71, 84, 84, 67, 71, 71, 71, 71, 84, 71, 84, 84, 84, 71, 71, 71, 84, 71, 65, 84, 65, 71, 65, 65, 65, 84, 84, 67, 84, 67, 71, 65, 65, 65, 84, 67, 67, 84, 84, 67, 84, 84, 65, 71, 71, 84, 67, 65, 65, 65, 65, 65, 65, 84, 65, 65, 65, 65, 65, 71, 84, 84, 71, 67, 71, 71, 84, 84, 71, 71, 71, 71], result: 0.002718
    ERROR ------------------------------------------------------------------------------------------
  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    OK, all I wanted was the GATK command line, but since it's in there... So, it looks like you're getting the dreaded "PairHMM Log Probability cannot be greater than 0" error, which I think is due to a known bug that we're working to fix. That said, we've only seen this bug with HaplotypeCaller. Are you saying you're also getting it with UnifiedGenotyper? If so, do you have an error log for that?

  • rfarrerrfarrer Cambridge MAMember

    Hi! No, I'm not getting that error with UG. But I am getting lines like this with UG:

    supercont_1.1 17505 . A . . . DP=2;MQ=60.00;MQ0=0 GT ./.

    supercont_1.1 17506 . T . . . DP=2;MQ=60.00;MQ0=0 GT ./.

    Should I assume these are also reference bases in addition to the following category of line (where genotype is given):

    supercont_1.1 1 . A . 73.23 . AN=2;DP=18;MQ=60.00;MQ0=0 GT:DP 0/0:17

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    No, GT ./. indicates that no call could be made because there was no useable coverage.

    In future, if you are experiencing separate issues with separate tools, please identify clearly which GATK command line resulted in which error message or problem. Or better yet, post them as separate questions. To be clear, posting a command line that runs a script, which itself may launch various disparate command lines, is not helpful for diagnosing issues.

  • rfarrerrfarrer Cambridge MAMember

    Thanks for the answer regarding coverage.

    Unfortunately, when it comes to the problems i'm having with Queue, I'm not sure how informative it would be to post the problem without including the script it is running, along with the command line that runs it? In any case, I didn't, and still do not, know where the problem is.

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    I'm sorry but I don't understand what you mean -- Queue seems to be working properly here. You got a GATK error message, which means that the problem lies with the GATK. Based on the command line of the GATK run, which is right there in the output log, you can see that your script is working fine (the command line makes perfect sense). The error message tells you that the HaplotypeCaller encountered an error while running PairHMM calculations. I told you this is a known bug that we are working on. What exactly is the problem you're having with Queue?

  • rfarrerrfarrer Cambridge MAMember

    OK. Well, GATK clearly gave an error message, but queue reported an error too (i.e. FunctionEdge - Error: 'java' ...).
    I wasn't sure if I had done something wrong, or either of the programs had, or any combination of the three.
    So, I will take it from your message that my qscript, queue, and my commands are fine, and the problem is with GATK.

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    Ah, sorry if that was unclear. Each individual job generated by Queue is a FunctionEdge in the job graph. So when you get a FunctionEdge error, it means that one of those jobs failed. Queue then gives you the details of which job failed, recapitulating the command line that was run, various parameters as applicable, and then the error/output of the program that was running in the failed job -- in this case the GATK HaplotypeCaller run.

    So yes, your script seems to work fine, and the issue is a bug in HaplotypeCaller that we are working to fix. For the record, I would recommend writing everything in a single QScript rather than generate/wrap different scripts within a shell script. But it may be a matter of personal preference because I don't use shell scripts much. Your approach is perfectly valid.

  • rfarrerrfarrer Cambridge MAMember

    In case it is of any use, I found that the problem went away when I removed ref-calls from the BQSR step. i.e. only calling refs on the final round of snp-calling.

Sign In or Register to comment.