QD in Join VCF vs its components in GVCFs does not match

Answers

• sum732 New YorkMember ✭

  February 2015

  Hello,
  To follow up
  Is this definition still holds true to calculate QD:
  QualByDepth (QD) calculation when used in large cohorts. Now, when the AD annotation is present for a given genotype then we only use its depth for QD if the variant depth > 1. Note that this only works in the gVCF pipeline for now.

  What if there is no AD value?

  For one of the samples I extracted all the "format" data and notice following:
  2614,072 GT:AB:AD:GQ:PL:SB
  2627,232 GT:AD:GQ:PL:SB
  123,63,8277 GT:DP:GQ:MIN_DP:PL
  63 GT:GQ:PL:SB
  So of all the total entries in a GVCF files majority does not have an "AD" value. How calculate QD then.

  Thanks

  0 · Share on Google+ Share on Twitter
• Sheila Broad InstituteMember, Broadie, Moderator admin

  February 2015

  @sum732
  Hi,

  Can you please confirm that you get the same results when using version 3.3?

  Thanks,
  Sheila

  0 · Share on Google+ Share on Twitter
• sum732 New YorkMember ✭

  February 2015

  Hello Sheila,
  Thanks for looking into this and responding.
  I have picked up this work done a few months ago. I don't think new version was used again, and ( if possible) newer version can be tested at a later stage.

  Do you think this is a bug? Is there any other explanation.
  To keep it simple, QD is based on AD Values only and it's calculated for each variant, across all samples and then normalized(the sum total is then used to divide the QUAL) is this correct.
  What happens when there is no AD?

  Best,

  0 · Share on Google+ Share on Twitter
• Geraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

  February 2015

  Older versions had some glitches with AD values so it could be a minor bug. But there may be some subtleties of the calculation that we need to clarify in the context of the reference confidence model.

  0 · Share on Google+ Share on Twitter
• sum732 New YorkMember ✭

  February 2015

  Hello Geraldine/Sheila,
  I had posted a comment yesterday and system said its "waiting approval". Request you to please release/approve/ that.

  For our analysis QD value (which a normalized val) is crucial. Without getting a complete understanding of it (how its calculated) cannot used it and can't move forward. Kindly elaborate on "subtleties of the calculation that we need to clarify in the context of the reference confidence model". What goes into it and how can I regenerate the values the match the report from GATK.
  Many Thanks.
  Best,

  0 · Share on Google+ Share on Twitter
• Sheila Broad InstituteMember, Broadie, Moderator admin

  February 2015

  @sum732

  Hi,

  I have asked two of our methods developers to clarify what exactly the calculation is, and I will get back to you when they respond.

  As for when there is no AD, you are looking at GVCF files. That is an intermediate file. The final vcf will have AD values for all sites.

  -Sheila

  0 · Share on Google+ Share on Twitter
• sum732 New YorkMember ✭

  February 2015

  Hi Shelia,

  Great. Will wait for their (->your) response.

  It sound bit strange that an AD for a variant if missing in intermediate files (GCVF), but will be present in the final output (joint VCF). So what value(s) are read to place a AD value for a sample if its missing in the individual and intermediate GVCF. I have also found that in some case there are no QD values for a variant call in joint VCFs (see below).

  Basically, I am retracing the steps to better understand the pipeline and for our analysis we need to identify the cause of difference and for that I need to be confident of all the values (calculations that goes into it). This is the reason I am finding these unexplained factors(? for a better word).

  e.g.(missing QD): AC=2;AF=0.071;AN=28;DP=590;FS=0.000;MQ=60.00;MQ0=0;Samples=Sample1;VQSLOD=18.23;VariantType=SNP;culprit=MQ
  e.g. (have):
  AC=2;AF=0.071;AN=28;DP=611;FS=0.000;MQ=60.00;MQ0=0;QD=30.59;Samples=Sample1;VQSLOD=19.91;VariantType=SNP;culprit=MQ

  Thanks for your understanding and patience.

  Best,

  0 · Share on Google+ Share on Twitter
• valentin Cambridge, MAMember, Dev ✭✭

  February 2015 Accepted Answer

  Hi @sum732,

  So as you noticed the QD is roughly speaking QUAL / Sum(Sum(AD)).

  Notice however that samples with no AD are ignored and so are samples that have AD but all the alternative alleles counts are zero. In other words QD does not take in consideration samples that do not show evidence for a concrete alternative allele. If no sample has a compliant AD then effectively QUAL is probably 0 (or very close) and so QD (probably not even present).

  However there are two important caveats here, two inaccuracy to the simple formula above.

  A.)

  First the Sum(Sum(AC)) is actually "normalized" using a weighted average length of the alternative alleles.

  So a more accurate formula for QD is

  QD = QUAL / (Sum(Sum(AD)) * alleleLengthNormalization(AC,ALT))

  Where alleleLengthNormalization(AC,ALT) = Max(1.0,weightedAverageAltAlleleLength(AC,ALT) / 3.0)

  Where weightedAverageAlleleLength(AC,ALT) = Sum(AC[i] * length(ALT[i]))/ Sum(AC).

  Where length(A) = is difference in number of bases between A allele and the reference and if this is zero the number of base mismatches between both sequences (thus for a SNP it would be 1, and for a simple insertion of 4 bases it would be 4).

  This normalization intent to reduce the QD for longer INDEL variants.

  B.)

  Actually QD has a soft-limit value of 35. If the QD calculated with the formula above is larger than 35, which happens amongst the examples that you provide; then this value is discarded and we draw a random QD from a normal distribution with mean 30 and std.dev. 3. That in turn means that several runs on the same data could potentially return different QDs although all close enough to 30.

  I'm afraid I can't answer the obvious question why 35/30/3 or why simply don't impose such a limit, perhaps another colleague of mine can add more to the discussion once is back from vacation.

  5 · Share on Google+ Share on Twitter
• valentin Cambridge, MAMember, Dev ✭✭

  February 2015

  @sum732 said:
  Hi Shelia,

  Great. Will wait for their (->your) response.

  It sound bit strange that an AD for a variant if missing in intermediate files (GCVF), but will be present in the final output (joint VCF). So what value(s) are read to place a AD value for a sample if its missing in the individual and intermediate GVCF.

  Looking at the examples GVCF that have around here there is no AD in hom-ref call blocks and there is AD when there is some significant evidence for a concrete alternative allele (i.e. anything except <NON-REF>). Do you have an example of a GVCF variant block with no AD?

  I have also found that in some case there are no QD values for a variant call in joint VCFs (see below).

  This is a bit more surprising to me... looking at those two very similar Info entries it seems strange that one would have it an the other not... it would be great if you could provide us with the full line that does not contain the QD and the steps to reproduce it (command and inputs focused on the problematic interval). This could be a bug.

  Basically, I am retracing the steps to better understand the pipeline and for our analysis we need to identify the cause of difference and for that I need to be confident of all the values (calculations that goes into it). This is the reason I am finding these unexplained factors(? for a better word).

  e.g.(missing QD): AC=2;AF=0.071;AN=28;DP=590;FS=0.000;MQ=60.00;MQ0=0;Samples=Sample1;VQSLOD=18.23;VariantType=SNP;culprit=MQ
  e.g. (have):
  AC=2;AF=0.071;AN=28;DP=611;FS=0.000;MQ=60.00;MQ0=0;QD=30.59;Samples=Sample1;VQSLOD=19.91;VariantType=SNP;culprit=MQ

  Thanks for your understanding and patience.

  Best,

  0 · Share on Google+ Share on Twitter
• sum732 New YorkMember ✭

  February 2015

  Hi Valentin.
  Thanks for the detailed response. I will try and post the entire line to help.

  0 · Share on Google+ Share on Twitter
• WimS Member ✭✭

  February 2016 edited February 2016

  Actually QD has a soft-limit value of 35. If the QD calculated with the formula above is larger than 35, which happens amongst the examples that you provide; then this value is discarded and we draw a random QD from a normal distribution with mean 30 and std.dev. 3. That in turn means that several runs on the same data could potentially return different QDs although all close enough to 30.

  I'm afraid I can't answer the obvious question why 35/30/3 or why simply don't impose such a limit, perhaps another colleague of mine can add more to the discussion once is back from vacation.

  I am curious what the answer is to the obvious 35/30/3 question. I also have sort of a custom situation where this formula might be of use. Is the colleague back who could answer this? Thank you.

  0 · Share on Google+ Share on Twitter

  Issue · Github

  February 2016 by Sheila

  Issue Number

  549

  State

  closed

  Last Updated

  February 2016

  Assignee

  Array

  Milestone

  Array

  Closed By

  vdauwera

• Geraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

  February 2016

  I'm not sure why the 35 cap; my guess is that this number was found empirically to work well for the purposes of variant filtering. The 30/3 part derives from 35 in that we want to introduce jitter into the value (to avoid having a huge pile of QDs of exactly 35) so that the annotation will behave well in VQSR, and again those specific numbers were probably found to work well to produce an adequate distribution of values.

  0 · Share on Google+ Share on Twitter
• WimS Member ✭✭

  February 2016

  @Geraldine_VdAuwera The link to VQSR seems to be probable. Thank you for the answer.

  0 · Share on Google+ Share on Twitter