search

dnanexus-rnd / GLnexus

Scalable gVCF merging and joint variant calling for population sequencing projects
Apache License 2.0
145 stars 37 forks source link

Bug in AF calculation when exome capture is used as --bed #168

Closed xunjieli closed 5 years ago

xunjieli commented 5 years ago

When exome capture regions is used as --bed, there are issues with AF calculation at the edges of these regions.

Here's one example where GLnexus will calculate the AF as 1.33333 when it should be 0.666667

HG002

19  13317490    .   GTTA    G,<*>   59  PASS    .   GT:GQ:DP:AD:VAF:PL  1/1:58:76:0,76,0:1,0:59,65,0,990,990,990

HG003

19  13317490    .   GTTA    G,<*>   69.7    PASS    .   GT:GQ:DP:AD:VAF:PL  0/1:66:76:34,42,0:0.552632,0:69,0,68,990,990,990

HG004

19  13317490    .   GTTA    G,<*>   70.7    PASS    .   GT:GQ:DP:AD:VAF:PL  0/1:69:78:37,41,0:0.525641,0:70,0,74,990,990,990

Bed file used (taken from agilent_sureselect_human_all_exon_v5_b37_targets.bed)

19  13317371    13317491
19  13317492    13317612
./glnexus_cli --config DeepVariant --bed chr19.bed HG00*.vcf.gz > merged.bcf

This merged.bcf will have AF=1.33333:

19  13317490    19_13317490_GTTA_G  GTTA    G   70  .   AF=1.33333;AQ=70    GT:DP:AD:GQ:PL:RNC  1/1:76:0,76:58:59,65,0:..   0/1:76:34,42:66:69,0,68:..  0/1:78:37,41:69:70,0,74:..

If repeat the above, but supply 19 0 10000000 as --bed, GLnexus is correct in the AF calculation. AF=0.666667:

19  13317490    19_13317490_GTTA_G  GTTA    G   70  .   AF=0.666667;AQ=70   GT:DP:AD:GQ:PL:RNC  1/1:76:0,76:58:59,65,0:..   0/1:76:34,42:66:69,0,68:..  0/1:78:37,41:69:70,0,74:..
mlin commented 5 years ago

:point_down: :man_facepalming: :point_down:

https://github.com/dnanexus-rnd/GLnexus/blob/master/src/discovery.cc#L71-L72

I think a possible fix will be to have merge_discovered_alleles recognize when the same deletion/MNP is being observed in multiple discovery ranges based on the discovered_allele_info::in_target field. Since the latter is optional, I need to double check that the latter field will be populated as expected in all cases.

xunjieli commented 5 years ago

Thanks, mlin@. The fix makes sense to me.

mlin commented 5 years ago

@xunjieli I implemented that fix on master, pls reopen if you still hit related issues. Expecting to cut a new point release next week after some other wip lands. Thanks!

xunjieli commented 5 years ago

Thanks @mlin. I've verified that this commit fixes the issue that I saw.