Overlapping positions

rfuentesrfuentes PhilippinesMember
edited February 2015 in Ask the GATK team

Hi

I'm wondering why some positions in VCF overlap. Can GATK skip emitting positions that are already part of an indel/concatenated ref?
We need to use EMIT_ALL_SITES but it's quite confusing if a position is represented multiple times especially if indel is present. I understand that there is always a duplicated position preceding an indel, but the positions after the indel should ideally be not overlapping the neighboring positions. Please see some examples below:

chr02 28792823 . C . 34.23 . AN=2;DP=12;MQ=57.31;MQ0=0 GT:DP 0/0:12
chr02 28792823 . CAA . 10000037.27 . AN=2;DP=12;MQ=57.31;MQ0=0 GT:AD:DP 0/0:12:12
chr02 28792824 . A . . . DP=12;MQ=57.31;MQ0=0 GT ./.
chr02 28792825 . A . 31.22 . AN=2;DP=12;MQ=57.31;MQ0=0 GT:DP 0/0:7

chr02 314879 . T . 100.23 . AN=2;DP=27;MQ=55.82;MQ0=0 GT:DP 0/0:27
chr02 314879 . TAGAG T,TAG 647.19 . AC=1,1;AF=0.500,0.500;AN=2;DP=27;FS=0.000;MLEAC=1,1;MLEAF=0.500,0.500;MQ=55.82;MQ0=0;QD=23.97;RPA=8,6,7;RU=AG;STR GT:AD:DP:GQ:PL 1/2:0,9,11:26:99:989,445,424,395,0,335
chr02 314880 . A . 43.23 . AN=2;DP=26;MQ=55.66;MQ0=0 GT:DP 0/0:5
chr02 314881 . G . 40.23 . AN=2;DP=25;MQ=55.47;MQ0=0 GT:DP 0/0:4
chr02 314882 . A . 73.23 . AN=2;DP=25;MQ=55.47;MQ0=0 GT:DP 0/0:16
chr02 314883 . G . 73.23 . AN=2;DP=25;MQ=55.47;MQ0=0 GT:DP 0/0:16

This next example even calls a SNP for position 10221400 even if deletion occurs in 10221399.

chr02 10221399 . C . 76.23 . AN=2;DP=17;MQ=58.74;MQ0=0 GT:DP 0/0:17
chr02 10221399 . CA CCTA,C 143.19 . AC=1,1;AF=0.500,0.500;AN=2;DP=17;FS=0.000;MLEAC=1,1;MLEAF=0.500,0.500;MQ=58.74;MQ0=0;QD=8.42 GT:AD:DP:GQ:PL 1/2:0,8,6:17:99:527,131,112,185,0,143
chr02 10221400 . A C,T 287.29 . AC=1,1;AF=0.500,0.500;AN=2;DP=17;Dels=0.29;FS=0.000;HaplotypeScore=6.7569;MLEAC=1,1;MLEAF=0.500,0.500;MQ=58.74;MQ0=0;QD=16.90 GT:AD:DP:GQ:PL 1/2:0,4,8:12:88:407,294,285,112,0,88
chr02 10221401 . A C 163.77 . AC=1;AF=0.500;AN=2;BaseQRankSum=0.854;DP=17;Dels=0.00;FS=1.913;HaplotypeScore=9.7562;MLEAC=1;MLEAF=0.500;MQ=58.74;MQ0=0;MQRankSum=1.558;QD=9.63;ReadPosRankSum=2.764 GT:AD:DP:GQ:PL 0/1:11,6:17:99:192,0,392
chr02 10221402 . A . 76.23 . AN=2;DP=17;MQ=58.74;MQ0=0 GT:DP 0/0:17

Thanks in advance for your help.
rfuentes

Answers

  • tommycarstensentommycarstensen United KingdomMember

    Hi @rfuentes . I'm a user like you. HaplotypeCaller merges SNPs and indels, but UnifiedGenotyper doesn't. Your variants at postions 10221399 (indel) and 10221400 (SNP) are different. Wouldn't leaving one of them out be like filtering your calls randomly?

    chr02 10221399 . C .
    chr02 10221399 . CA CCTA,C
    chr02 10221400 . A C,T
    
  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie

    To clarify further, the reason why UG does this is because it uses separate modeling algorithms to call SNPs and indels. The program doesn't actually "know" that it is calling a SNP and an indel at the same position. Obviously that's not optimal, which is why in HaplotypeCaller we made the program smarter and capable of evaluating SNPs and indels at the same time. So HC does not emit them separately. You should really consider using HC instead of UG, it's much better in many ways.

  • rfuentesrfuentes PhilippinesMember

    Thank you @tommycarstensen and @Geraldine_VdAuwera! I will inform my colleague about this. I'm not sure why they chose to use UG in the pipeline instead of HC. Do you think the VCF is unreliable because of these mixed cases or can you suggest any method so we can select the right call?

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie

    UG was the standard fo long time, and some people have not yet updated their pipeline, it's understandable. We do not provide recommendations for choosing in mixed cases since it can depend on project objectives. But perhaps @tommycarstensen can share his experiences... :)

  • tommycarstensentommycarstensen United KingdomMember

    If this is low (or high) coverage data and you have at least 100 samples, then you should see reasonably good sensitivity and specificity with UG. I would find the calls reliable. Maybe you have some SNP array data or high coverage data for the same samples, which you can compare against?

Sign In or Register to comment.