GenotypeGVCFs not outputting some variants it should with -stand-call-conf
open | Created 2019-03-13 | Last updated 2019-04-16| Posted by tfenne | See in Github


VariantCalling bug


Bug Report

Affected tool(s) or class(es)

GenotypeGVCFs 4.0.0.12

Affected version(s)

  • Latest public release version [version?]
  • Latest master branch as of [date of test?]

Description

I've run into a weird case where GenotypeGVCFs is doing something unexpected. I have a gVCF with the following entry in it:

chr11   6637739 .       ATTTTT  A,AT,ATT,ATTT,ATTTT,ATTTTTT,<NON_REF>   565.73  .       BaseQRankSum=-0.014;ClippingRankSum=0.508;DP=94;ExcessHet=3.0103;MLEAC=0,0,0,1,0,0,0;MLEAF=0,0,0,0.5,0,0,0;MQRankSum=0;RAW_MQandDP=338400,94;REF_BASES=GCCGGCCTGGATTTTTTTTTT;ReadPosRankSum=-0.812      GT:AD:DP:F1R2:F2R1:GQ:PL:SB     0/4:9,3,3,11,15,8,3,0:52:8,2,2,8,12,6,3,0:1,1,1,3,3,2,0,0:56:603,504,1526,335,1171,1118,56,661,640,608,0,362,313,183,335,336,500,389,187,171,527,655,864,622,277,169,466,1026,597,1101,953,645,465,625,861,1133:8,1,33,10

It's a messy site for sure, an indel in a long homopolymer-T, but I think that's a separate issue. If I run the following on that gVCF:

gatk GenotypeGVCFs \
  -R hg19.fa -V test.g.vcf -O test.vcf \
  -A ClippingRankSumTest -A Coverage -A ExcessHet -A FisherStrand \
  -A MappingQualityRankSumTest -A OxoGReadCounts -A QualByDepth -A ReadPosRankSumTest \
  -A ReferenceBases -A RMSMappingQuality -A StrandOddsRatio -A TandemRepeat \
  -L chr11:6637730-6637750 \
  -stand-call-conf 18.0 \

then I get the following output to the VCF just like I'd expect:

chr11   6637739 .       ATT     A       565.73  .       AC=1;AF=0.500;AN=2;BaseQRankSum=-1.400e-02;ClippingRankSum=0.508;DP=94;ExcessHet=3.0103;FS=1.779;MLEAC=1;MLEAF=0.500;MQ=60.00;MQRankSum=0.00;QD=23.57;REF_BASES=GCCGGCCTGGATTTTTTTTTT;RPA=15,13;RU=T;ReadPosRankSum=-8.120e-01;SOR=0.386;STR    GT:AD:DP:F1R2:F2R1:GQ:PL        0/1:9,15:52:8,2,2,8,12,6,3,0:1,1,1,3,3,2,0,0:99:603,0,335

QUAL is unchanged since I'm genotyping a single-sample gVCF. However, if I raise my -stand-call-conf threshold to 19.0, GenotypeGVCFs no longer outputs any variants. 565.73 >> 19.0, so I'm confused as to why that variant is no longer emitted.

Steps to reproduce

Run GenotypeGVCFs on the above example with -stand-call-conf values ranging up to ~550.

Expected behavior

The variant should be emitted into the VCF.

Actual behavior

The variant is not emitted if -stand-call-conf is set to 19 or higher.


Return to top