Make CombineGVCFs and GenomicsDBImport properly compute AD over overlapping deletions
open | Created 2018-09-05 | Last updated 2019-02-05| Posted by cwhelan | See in Github


HaplotypeCaller bug


When #4963 is merged, HaplotypeCaller will produce more spanning deletions, and in some instances will produce spanning deletion alleles for deletions that are overlapped by other deletions. In those instances, CombineGVCFs and GenomicsDBImport do not properly compute the genotype FORMAT AD and DP tags for sites overlapped by the deletions.

For example, for the input GVCF lines

20	10068158	.	GTGTATATATATA	G,<NON_REF>	66.73	.	BaseQRankSum=-0.652;ClippingRankSum=0.000;DP=29;ExcessHet=3.0103;MLEAC=1,0;MLEAF=0.500,0.00;MQRankSum=0.328;RAW_MQ=93364.00;ReadPosRankSum=-0.253	GT:AD:DP:GQ:PL:SB	0/1:3,4,0:7:57:104,0,57,114,69,183:0,3,2,2
20	10068160	.	GTATATATATATGTA	G,*,<NON_REF>	697.73	.	DP=28;ExcessHet=3.0103;MLEAC=1,1,0;MLEAF=0.500,0.500,0.00;RAW_MQ=87005.00	GT:AD:DP:GQ:PL:SB	1/2:0,2,4,0:6:53:735,162,131,102,0,53,507,174,108,472:0,0,2,4

Combine GVCFs run as follows:

./gatk CombineGVCFs -V src/test/resources/org/broadinstitute/hellbender/tools/haplotypecaller/expected.testGVCFMode.gatk4.g.vcf -O test_gdb_import_combine.g.vcf -R src/test/resources/large/human_g1k_v37.20.21.fasta

Returns the following output:

20      10068158        .       GTGTATATATATA   G,<NON_REF>     .       .       BaseQRankSum=-6.520e-01;ClippingRankSum=0.00;DP=29;ExcessHet=3.01;MQRankSum=0.328;RAW_MQ=93364.00;ReadPosRankSum=-2.530e-01     GT:AD:DP:GQ:PL:SB       ./.:3,4,0:7:57:104,0,57,114,69,183:0,3,2,2
20      10068159        .       T       *,<NON_REF>     .       .       DP=29   GT:AD:DP:GQ:PL:SB       ./.:3,4,0:7:57:104,0,57,114,69,183:0,3,2,2
20      10068160        .       GTATATATATATGTA G,*,<NON_REF>   .       .       DP=28;ExcessHet=3.01;RAW_MQ=87005.00    GT:AD:DP:GQ:PL:SB       ./.:0,2,4,0:6:53:735,162,131,102,0,53,507,174,108,472:0,0,2,4
20      10068161        .       T       *,<NON_REF>     .       .       DP=28   GT:AD:DP:GQ:PL:SB       ./.:0,4,0:6:53:735,102,53,507,108,472:0,0,2,4
20      10068162        .       A       *,<NON_REF>     .       .       DP=28   GT:AD:DP:GQ:PL:SB       ./.:0,4,0:6:53:735,102,53,507,108,472:0,0,2,4
20      10068163        .       T       *,<NON_REF>     .       .       DP=28   GT:AD:DP:GQ:PL:SB       ./.:0,4,0:6:53:735,102,53,507,108,472:0,0,2,4
20      10068164        .       A       *,<NON_REF>     .       .       DP=28   GT:AD:DP:GQ:PL:SB       ./.:0,4,0:6:53:735,102,53,507,108,472:0,0,2,4
20      10068165        .       T       *,<NON_REF>     .       .       DP=28   GT:AD:DP:GQ:PL:SB       ./.:0,4,0:6:53:735,102,53,507,108,472:0,0,2,4
20      10068166        .       A       *,<NON_REF>     .       .       DP=28   GT:AD:DP:GQ:PL:SB       ./.:0,4,0:6:53:735,102,53,507,108,472:0,0,2,4
20      10068167        .       T       *,<NON_REF>     .       .       DP=28   GT:AD:DP:GQ:PL:SB       ./.:0,4,0:6:53:735,102,53,507,108,472:0,0,2,4
20      10068168        .       A       *,<NON_REF>     .       .       DP=28   GT:AD:DP:GQ:PL:SB       ./.:0,4,0:6:53:735,102,53,507,108,472:0,0,2,4
20      10068169        .       T       *,<NON_REF>     .       .       DP=28   GT:AD:DP:GQ:PL:SB       ./.:0,4,0:6:53:735,102,53,507,108,472:0,0,2,4
20      10068170        .       A       *,<NON_REF>     .       .       DP=28   GT:AD:DP:GQ:PL:SB       ./.:0,4,0:6:53:735,102,53,507,108,472:0,0,2,4
20      10068171        .       T       *,<NON_REF>     .       .       DP=28   GT:AD:DP:GQ:PL:SB       ./.:0,4,0:6:53:735,102,53,507,108,472:0,0,2,4
20      10068172        .       G       *,<NON_REF>     .       .       DP=28   GT:AD:DP:GQ:PL:SB       ./.:0,4,0:6:53:735,102,53,507,108,472:0,0,2,4
20      10068173        .       T       *,<NON_REF>     .       .       DP=28   GT:AD:DP:GQ:PL:SB       ./.:0,4,0:6:53:735,102,53,507,108,472:0,0,2,4
20      10068174        .       A       *,<NON_REF>     .       .       DP=28   GT:AD:DP:GQ:PL:SB       ./.:0,4,0:6:53:735,102,53,507,108,472:0,0,2,4

GenomicsDBImport run like this:

./gatk GenomicsDBImport -R src/test/resources/large/human_g1k_v37.20.21.fasta -L 20 -V src/test/resources/org/broadinstitute/hellbender/tools/haplotypecaller/expected.testGVCFMode.gatk4.g.vcf -genomicsdb-workspace-path spanDelWorkspace
./gatk SelectVariants -V gendb://spanDelWorkspace -R src/test/resources/large/human_g1k_v37.20.21.fasta -O test.g.vcf -L 20

Returns the following output:

20      10068158        .       GTGTATATATATA   G,<NON_REF>     .       .       BaseQRankSum=-6.520e-01;ClippingRankSum=0.00;DP=29;ExcessHet=3.01;MQRankSum=0.328;RAW_MQ=93364.00;ReadPosRankSum=-2.530e-0
1     GT:AD:DP:GQ:PL:SB       ./.:3,4,0:7:57:104,0,57,114,69,183:0,3,2,2
20      10068159        .       T       *,<NON_REF>     .       .       DP=29   GT:AD:DP:GQ:PL:SB       ./.:3,4,0:7:57:104,0,57,114,69,183:0,3,2,2
20      10068160        .       GTATATATATATGTA G,*,<NON_REF>   .       .       DP=28;ExcessHet=3.01;RAW_MQ=87005.00    GT:AD:DP:GQ:PL:SB       ./.:0,2,4,0:6:53:735,162,131,102,0,53,507,174,108,472:0,0,
2,4
20      10068161        .       T       *,<NON_REF>     .       .       DP=28   GT:AD:DP:GQ:PL:SB       ./.:0,2,0:6:53:735,162,131,507,174,472:0,0,2,4
20      10068162        .       A       *,<NON_REF>     .       .       DP=28   GT:AD:DP:GQ:PL:SB       ./.:0,2,0:6:53:735,162,131,507,174,472:0,0,2,4
20      10068163        .       T       *,<NON_REF>     .       .       DP=28   GT:AD:DP:GQ:PL:SB       ./.:0,2,0:6:53:735,162,131,507,174,472:0,0,2,4
20      10068164        .       A       *,<NON_REF>     .       .       DP=28   GT:AD:DP:GQ:PL:SB       ./.:0,2,0:6:53:735,162,131,507,174,472:0,0,2,4
20      10068165        .       T       *,<NON_REF>     .       .       DP=28   GT:AD:DP:GQ:PL:SB       ./.:0,2,0:6:53:735,162,131,507,174,472:0,0,2,4
20      10068166        .       A       *,<NON_REF>     .       .       DP=28   GT:AD:DP:GQ:PL:SB       ./.:0,2,0:6:53:735,162,131,507,174,472:0,0,2,4
20      10068167        .       T       *,<NON_REF>     .       .       DP=28   GT:AD:DP:GQ:PL:SB       ./.:0,2,0:6:53:735,162,131,507,174,472:0,0,2,4
20      10068168        .       A       *,<NON_REF>     .       .       DP=28   GT:AD:DP:GQ:PL:SB       ./.:0,2,0:6:53:735,162,131,507,174,472:0,0,2,4
20      10068169        .       T       *,<NON_REF>     .       .       DP=28   GT:AD:DP:GQ:PL:SB       ./.:0,2,0:6:53:735,162,131,507,174,472:0,0,2,4
20      10068170        .       A       *,<NON_REF>     .       .       DP=28   GT:AD:DP:GQ:PL:SB       ./.:0,2,0:6:53:735,162,131,507,174,472:0,0,2,4
20      10068171        .       T       *,<NON_REF>     .       .       DP=28   GT:AD:DP:GQ:PL:SB       ./.:0,2,0:6:53:735,162,131,507,174,472:0,0,2,4
20      10068172        .       G       *,<NON_REF>     .       .       DP=28   GT:AD:DP:GQ:PL:SB       ./.:0,2,0:6:53:735,162,131,507,174,472:0,0,2,4
20      10068173        .       T       *,<NON_REF>     .       .       DP=28   GT:AD:DP:GQ:PL:SB       ./.:0,2,0:6:53:735,162,131,507,174,472:0,0,2,4
20      10068174        .       A       *,<NON_REF>     .       .       DP=28   GT:AD:DP:GQ:PL:SB       ./.:0,2,0:6:53:735,162,131,507,174,472:0,0,2,4
20      10068175        .       T       <NON_REF>       .       .       .       GT:DP:GQ:MIN_DP:PL      ./.:20:17:20:0,17,729

In this case there are two overlapping deletions, one spanning from 10068158-10068170 and the other spanning from 10068160-10068174.

AD and DP should be as follows for each site:

10068158 3,4,0:7
10068159 3,4,0:7
10068160 0,2,4,0:6
10068161 0,6,0:6
10068162 0,6,0:6
10068163 0,6,0:6
10068164 0,6,0:6
10068165 0,6,0:6
10068166 0,6,0:6
10068167 0,6,0:6
10068168 0,6,0:6
10068179 0,6,0:6
10068170 0,6,0:6
10068171 4,2,0:6
10068172 4,2,0:6
10068173 4,2,0:6
10068174 4,2,0:6

For the final set of sites that are overlapped by the second deletion but not the first, the reads that supported the first deletion should return to being counted as reference-supporting reads.


Return to top