Incorrect GQ Score = 0
closed | Created 2018-11-24 | Last updated 2019-08-02| Posted by jjfarrell | See in Github


HaplotypeCaller bug


Bug Report

Affected tool(s) or class(es)

HaplotypeCaller

Affected version(s)

GATK 3.7 and GATK 4.0.11.0

Description

GATK Pipeline (HaplotypeCaller->gvcf-->importdb->GenotypeGVCF)is calling sites with a GQ=0. But these sites often have plenty of coverage and no obvious reason for such a low GQ score. Often the GQ should be 99 as the DP >40. This seems to be primarily an issue with homozygous reference calls.

The GT is accurate for the high DP sites but the inaccurate GQ is problematic for any genotype level qc on the pVCF. If the site is recoded from 0/0 to './.' for GQ <20, the result is higher missing rate due to the inaccurate GQ=0.

Directly calling the VCF with HaplotypeCaller without the gVCF intermediate gVCF file calculates the correct GQ score. Freebayes also calculates a correct GQ on these samples.
rs429358_gq_dp.pdf

Steps to reproduce

I am seeing this bug for 57 samples of 5000 crams at snp rs429358 but I would expect it is not unique to this site.

Select two crams with a Passed site with:
cram 1. Call with GT='0/0, GQ=0 and DP >40.
cram 2. Call with GT='0/1' or '1/1' and DP>20.

Create vcf with two approaches:

Pipeline 1. HaplotypeCaller-->vcf

module load gatk/4.0.11.0
gatk HaplotypeCaller -R /restricted/projectnb/casa/ref/GRCh38_full_analysis_set_plus_decoy_hla.fa
-I gq0_cram.list
-L chr19:44907684-44909822
--use-new-qual-calculator
-O good.vcf.gz

Good GQ scores were also estimated with Freebayes on these samples also.

Pipeline 2 HaplotypeCaller --> bvcf--->ImportVCF-->GenotypeVCF-->VCF with 2 samples

gatk HaplotypeCaller -R /restricted/projectnb/casa/ref/GRCh38_full_analysis_set_plus_decoy_hla.fa
-I $sample.cram
--use-new-qual-calculator
-L chr19:44907684-44909822
-ERC GVCF
-O bad.g.vcf.gz

Followed by import and GenotypeVCF.

Expected behavior

Pipeline 2 should generate accurate GQ scores that match the GQ in the HaplotypeCaller vcf output of pipeline 1. Instead GQ=0.

This is the output for the 57 GQ=0 samples with pipeline 1 which is accurate.

AC=7;AF=0.061;AN=114;BaseQRankSum=-6.147;DP=1846;ExcessHet=3.8592;FS=0.000;InbreedingCoeff=-0.0640;MLEAC=6;MLEAF=0.053;MQ=60.00;MQRankSum=0.000;QD=4.52;ReadPosRankSum=-0.781;SOR=2.833
GT:AD:DP:GQ:PL
0/0:37,0:37:99:0,111,1236
0/0:40,0:40:99:0,120,1357
0/0:34,0:34:99:0,102,1161
0/0:49,0:49:99:0,147,1673
0/0:33,0:33:99:0,99,1036
0/0:48,0:48:99:0,144,1728
0/0:42,0:42:99:0,126,1410
0/0:37,0:37:99:0,111,1215
0/0:39,0:39:99:0,117,1311
0/0:42,0:42:99:0,126,1419
0/0:53,0:53:99:0,159,1744
0/0:45,0:45:99:0,135,1529
0/0:44,0:44:99:0,132,1419
0/0:38,0:38:99:0,114,1299
0/0:37,0:37:99:0,111,1205
0/0:34,0:34:99:0,102,1151
0/0:57,0:57:99:0,171,1826
0/0:27,1:28:49:0,49,904
0/0:41,0:41:99:0,123,1364
0/0:28,0:28:84:0,84,933
0/0:36,0:36:99:0,108,1171
0/0:29,0:29:87:0,87,987
0/0:31,0:31:93:0,93,997
0/0:37,0:37:99:0,111,1266
0/0:28,2:30:70:0,70,914
0/0:36,0:36:99:0,108,1230
0/0:49,0:49:99:0,147,1613
0/0:38,1:39:82:0,82,1231
0/0:26,0:26:78:0,78,869
0/0:22,0:22:66:0,66,734
0/0:28,1:29:52:0,52,878
0/0:36,0:36:99:0,108,1221
0/0:39,0:39:99:0,117,1315
0/0:31,1:32:86:0,86,1008
0/0:38,0:38:99:0,114,1286
0/0:37,0:37:99:0,111,1266
0/0:35,0:35:99:0,105,1131
0/0:32,0:32:96:0,96,1071
0/0:46,0:46:99:0,138,1508
0/0:27,0:27:81:0,81,918
0/0:19,0:19:57:0,57,620
0/0:45,0:45:99:0,135,1474
0/0:50,0:50:99:0,150,1687
0/0:29,0:29:87:0,87,968
0/0:30,0:30:90:0,90,992
0/0:55,1:56:99:0,158,1830
0/1:8,2:10:19:19,0,248
0/1:5,1:6:17:17,0,145
0/1:4,1:5:10:10,0,124
0/1:5,1:6:17:17,0,155
0/1:3,2:5:34:34,0,79
0/1:5,3:8:45:45,0,150
0/1:2,1:3:26:26,0,61
0/0:13,0:13:39:0,39,431
0/0:17,0:17:51:0,51,571
0/0:28,0:28:84:0,84,993
0/0:10,0:10:30:0,30,328

Freebayes has been run on these 57 samples and also get '0/0' but with GQ in the 140-160 range for most samples and are in line with the results with the HaplotypeCaller VCF direct results (see below).

Actual behavior

Pipeline 2 is incorrectly setting GQ=0.

This is the output for the 57 GQ=0 samples+ 1 extra sample GQ=99 with pipeline 2. The extra sample is needed to produce a variant record otherwise all the records would be homoz refer with GQ=0.

AC=1;AF=8.621e-03;AN=116;BaseQRankSum=-8.310e-01;DP=2213;ExcessHet=41.0061;FS=1.957;InbreedingCoeff=-0.3410;MLEAC=10;MLEAF=0.086;MQ=60.00;MQRankSum=0.00;QD=12.17;ReadPosRankSum=0.616;SOR=1.080
GT:AD:DP:GQ:PL
0/1:9,8:17:99:227,0,272
0/0:41,0:41:0:0,0,1097
0/0:51,0:51:0:0,0,1216
0/0:61,0:61:0:0,0,1373
0/0:54,0:54:0:0,0,962
0/0:49,0:49:0:0,0,1156
0/0:53,0:53:0:0,0,729
0/0:44,0:44:0:0,0,1161
0/0:38,0:38:0:0,0,963
0/0:68,0:68:0:0,0,1518
0/0:33,0:33:0:0,0,841
0/0:54,0:54:0:0,0,687
0/0:44,0:44:0:0,0,1003
0/0:33,0:33:0:0,0,709
0/0:54,0:54:0:0,0,580
0/0:31,0:31:0:0,0,790
0/0:36,0:36:0:0,0,843
0/0:49,0:49:0:0,0,978
0/0:34,0:34:0:0,0,669
0/0:39,0:39:0:0,0,898
0/0:60,0:60:0:0,0,1270
0/0:48,0:48:0:0,0,908
0/0:30,0:30:0:0,0,780
0/0:44,0:44:0:0,0,778
0/0:24,0:24:0:0,0,531
0/0:40,0:40:0:0,0,714
0/0:50,0:50:0:0,0,794
0/0:44,0:44:0:0,0,953
0/0:38,0:38:0:0,0,953
0/0:45,0:45:0:0,0,1081
0/0:40,0:40:0:0,0,967
0/0:49,0:49:0:0,0,604
0/0:38,0:38:0:0,0,675
0/0:55,0:55:0:0,0,1111
0/0:65,0:65:0:0,0,1056
0/0:32,0:32:0:0,0,734
0/0:21,0:21:0:0,0,551
0/0:51,0:51:0:0,0,1335
0/0:56,0:56:0:0,0,1549
0/0:40,0:40:0:0,0,829
0/0:35,0:35:0:0,0,848
0/0:63,0:63:0:0,0,1409
0/0:10,0:10:0:0,0,229
0/0:6,0:6:0:0,0,129
0/0:5,0:5:0:0,0,114
0/0:38,0:38:0:0,0,893
0/0:6,0:6:0:0,0,139
0/0:5,0:5:0:0,0,45
0/0:9,0:9:0:0,0,75
0/0:3,0:3:0:0,0,35
0/0:16,0:16:0:0,0,199
0/0:22,0:22:0:0,0,526
0/0:32,0:32:0:0,0,924
0/0:13,0:13:0:0,0,313
0/0:55,0:55:0:0,0,1394
0/0:50,0:50:0:0,0,1294
0/0:45,0:45:0:0,0,863
0/0:47,0:47:0:0,0,1092

The problem appears to be occurring in the generation of the gVCF. The GQ=0 is being output for that region in the gVCF which is propagated to the VCF output. Here are some of the gVCF records for this site.

chr19 44908684 . T <NON_REF> . . END=44908688 GT:DP:GQ:MIN_DP:PL 0/0:44:0:41:0,0,1097
chr19 44908684 . T <NON_REF> . . END=44908688 GT:DP:GQ:MIN_DP:PL 0/0:47:0:44:0,0,1003
chr19 44908684 . T <NON_REF> . . END=44908688 GT:DP:GQ:MIN_DP:PL 0/0:44:0:44:0,0,778
chr19 44908684 . T <NON_REF> . . END=44908688 GT:DP:GQ:MIN_DP:PL 0/0:66:0:65:0,0,1056
chr19 44908684 . T <NON_REF> . . END=44908688 GT:DP:GQ:MIN_DP:PL 0/0:39:0:38:0,0,893
chr19 44908684 . T <NON_REF> . . END=44908688 GT:DP:GQ:MIN_DP:PL 0/0:55:0:55:0,0,1394
chr19 44908684 . T <NON_REF> . . END=44908688 GT:DP:GQ:MIN_DP:PL 0/0:51:0:50:0,0,1294
chr19 44908684 . T <NON_REF> . . END=44908688 GT:DP:GQ:MIN_DP:PL 0/0:46:0:45:0,0,863
chr19 44908684 . T <NON_REF> . . END=44908688 GT:DP:GQ:MIN_DP:PL 0/0:48:0:47:0,0,1092
chr19 44908684 . T <NON_REF> . . END=44908688 GT:DP:GQ:MIN_DP:PL 0/0:53:0:51:0,0,1216


Return to top