Holiday Notice:
The Frontline Support team will be slow to respond December 17-18 due to an institute-wide retreat and offline December 22- January 1, while the institute is closed. Thank you for your patience during these next few weeks. Happy Holidays!

Possible bug in SelectVariants tool

Dear GATK experts,

I have done variant calling on 384 potato samples following, mostly, GATK best practices and have applied hard filters to select SNPs for further usage. However, I am noticing that '--max-nocall-fraction', '--max-nocall-number' and '--max-fraction-filtered-genotypes' arguments for 'SelectVariants' are not working properly. I have tried with various cutoff settings and every time I am observing SNPs with a much larger number of genotypes (~246 out of 384) with 'no call' than the set thresholds. I have searched the forum first but couldn't find any relevant threads. I am using the latest GATK version (4.0.7.0). I am attaching three example sets of (1) log files (2) subset vcf files and (3) vcf index file for the three main vcfs. I would appreciate if you could provide any feedback on this issue and/or if this behaviour has been observed by some other users also.

Regards,
Sanjeev

Issue · Github
by Sheila

Issue Number
3147
State
open
Last Updated
Assignee
Array

Best Answer

Answers

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @sanjeevksh
    Hi Sanjeev,

    I will take a look and get back to you.

    -Sheila

  • sanjeevkshsanjeevksh Member

    @Sheila
    Thanks for taking up this request,
    Sanjeev

  • sanjeevkshsanjeevksh Member

    @Sheila
    Is there any update on the issue I mentioned?
    Thanks,
    Sanjeev

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @sanjeevksh
    Hi Sanjeev,

    No, I have not had a chance to look at this yet. I will try to take a look in the next few weeks and get back to you.

    -Sheila

  • sanjeevkshsanjeevksh Member

    @Sheila
    Thanks, that would be really great,
    Regards,
    Sanjeev

  • sanjeevkshsanjeevksh Member

    @Sheila
    Hi Sheila,
    I am wondering if you have had any chance to look into the issue I mentioned, would highly appreciate some early help.
    Regards,
    Sanjeev

  • bhanuGandhambhanuGandham Member, Administrator, Broadie, Moderator admin

    Hi @sanjeevksh

    I will be taking over for sheila.
    I have asked the developers for an update on this and I will get back to you soon.
    Sorry for the delay.

    Regards
    Bhanu

  • sanjeevkshsanjeevksh Member

    Hi @bhanuGandham
    Many thanks for this update and I really appreciate your efforts to expedite the progress with this request.
    Regards,
    Sanjeev

  • bhanuGandhambhanuGandham Member, Administrator, Broadie, Moderator admin

    Hi @sanjeevksh

    Would you please send me the exact command you are using and a few chr locations you are seeing this discrepancy. thank you.

    Regards
    Bhanu

  • sanjeevkshsanjeevksh Member
    edited November 13

    Hi @bhanuGandham
    Sorry I saw your message today only.
    Here is the command I used:

    gatk SelectVariants \
    -R ~/reference_DM_PM4.03_G3_all_bowtie2/DM_v4.03_G3_all.fasta \
    -V $1.AllChr.final.fltd.variants.QD10.SOR5.MQ20.DP20t1000.GQ20.vcf.gz \
    --exclude-filtered \
    --select-type-to-include SNP \
    --restrict-alleles-to BIALLELIC \
    --set-filtered-gt-to-nocall \
    --max-nocall-number 42 \
    --selectExpressions 'vc.getGenotype("GPP116880").isHomRef()' \
    --exclude-non-variants \
    --add-output-vcf-command-line \
    -O $1.AllChr.final.fltd.slctd.variants.QD10.SOR5.MQ20.DP20t1000.GQ20.noCallGeno42.vcf.gz

    There is no specific chromosome location as such where this issue occurs more predominantly.

    The main issue is that despite setting cutoff settings for removing SNPs with 'no call', SNPs with a much larger number of missing calls than the set threshold are present in the SelectVariant processed vcf. For example, with a setting of '--max-nocall-number 42' SNPs with missing calls for ~246 genotypes (out of 384) are still present.

    I have tried all three arguments ('--max-nocall-fraction', '--max-nocall-number' and '--max-fraction-filtered-genotypes') to filter on 'missing calls' threshold but none of them seem to work.

    My vcf has different ploidy (diploid, triploid, tetraploid and hexaploid) samples, do you think this could be causing the issue?

    Regards,
    Sanjeev

    Post edited by sanjeevksh on
  • bhanuGandhambhanuGandham Member, Administrator, Broadie, Moderator admin

    Hi @sanjeevksh

    What i mean is would you please post a few records from your vcf where you see that nocall present in more than 42 samples? This will help us debug.

    Regards
    Bhanu

  • sanjeevkshsanjeevksh Member

    Hi @bhanuGandham

    Actually I noticed all SNPs in the subset file I sent to you have nocall for more than 42 samples ranging from 55 to 327.

    Top twenty locations with highest nocall samples are as below:

    CHROM POS ID REF ALT QUAL FILTER #missingGeno

    ST4.03ch01 13752067 . A T 5976.83 PASS 327
    ST4.03ch01 13752090 . G A 1159.59 PASS 323
    ST4.03ch01 13752116 . G A 13239.34 PASS 322
    ST4.03ch01 13752117 . G A 2827.22 PASS 322
    ST4.03ch01 13752068 . G A 20084.84 PASS 321
    ST4.03ch01 13752065 . G C 745.65 PASS 320
    ST4.03ch01 13752108 . C T 421.76 PASS 320
    ST4.03ch01 13752119 . G A 1780 PASS 320
    ST4.03ch01 16249840 . C A 3528.45 PASS 294
    ST4.03ch01 13752149 . T G 4244.97 PASS 261
    ST4.03ch01 13752156 . C T 226.02 PASS 259
    ST4.03ch01 18393338 . T C 13227.22 PASS 246
    ST4.03ch01 18393388 . T C 78453.73 PASS 243
    ST4.03ch01 18393370 . G A 3782.04 PASS 242
    ST4.03ch01 18393424 . A G 598.93 PASS 241
    ST4.03ch01 18393425 . T A 3422.12 PASS 241
    ST4.03ch01 18393413 . T C 853.96 PASS 239
    ST4.03ch01 18393335 . T C 1002.7 PASS 238
    ST4.03ch01 18393347 . G A 10907.31 PASS 238
    ST4.03ch01 18393359 . C T 800.8 PASS 238

    Twenty locations with lowest nocall samples are as below:

    CHROM POS ID REF ALT QUAL FILTER #missingGeno

    ST4.03ch01 6196053 . T C 628.52 PASS 76
    ST4.03ch01 6196055 . T C 603.25 PASS 76
    ST4.03ch01 6196082 . C T 4199.82 PASS 76
    ST4.03ch01 6196077 . T C 1473.17 PASS 75
    ST4.03ch01 6196079 . C T 2130.93 PASS 75
    ST4.03ch01 6196083 . T C 1166.52 PASS 75
    ST4.03ch01 6295515 . T C 259.21 PASS 75
    ST4.03ch01 6196061 . G A 3101.23 PASS 74
    ST4.03ch01 6196087 . G T 594.25 PASS 74
    ST4.03ch01 6196107 . A G 6124.66 PASS 74
    ST4.03ch01 6196109 . T C 513.26 PASS 74
    ST4.03ch01 6196132 . C T 544.26 PASS 74
    ST4.03ch01 6196141 . G C 927.26 PASS 74
    ST4.03ch01 6196143 . C T 2187.14 PASS 74
    ST4.03ch01 309589 . G T 105695.55 PASS 73
    ST4.03ch01 6302470 . C T 745.51 PASS 65
    ST4.03ch01 18393200 . A C 3384.28 PASS 64
    ST4.03ch01 18393213 . C T 548.93 PASS 62
    ST4.03ch01 309580 . T G 929.26 PASS 57
    ST4.03ch01 4018671 . A G 39948.25 PASS 55

    Please let me know if you need any further details,

    Regards,
    Sanjeev

  • bhanuGandhambhanuGandham Member, Administrator, Broadie, Moderator admin

    Hi @sanjeevksh

    we are looking into this issue and will get back to you soon.

    Regards
    Bhanu

  • bhanuGandhambhanuGandham Member, Administrator, Broadie, Moderator admin

    Hi @sanjeevksh

    This is the result of multiple options being invoked and not activating in the order we expect.
    I suspect that filtered samples will be set to "no call" and then if the total of true "no call" samples plus the filtered is >42, then the whole site will be excluded. In the code however, --max-nocall-number is evaluated before filtered samples are set to "no call". So I suspect that what's happening is:
    1) None of these sites have > 42 "no call" samples, so they pass that filter
    2) However some of them have filtered genotypes, which are subsequently converted to "no call"
    3) The final vcf appears to have more than 42 total "no call" samples in many sites.

    If this is correct, the fix would be to
    1) don't use --set-filtered-gt-to-nocall (at least at first, for diagnostic purposes)
    2) use --max-nocall-number 21 --max-filtered-genotypes 21
    In which case we would expect to see the offending sites filtered out.

    Please let me know how this suggestion works out.

    Thank you.

    Regards
    Bhanu

  • sanjeevkshsanjeevksh Member
    edited November 22

    Hi @bhanuGandham

    Thanks for your feedback.

    Few points/clarifications as detailed below:

    (1) Do you mean the order of activating different options is not correct in the command I submitted or the way this tool is set up to process the different options supplied by the user? In the command I submitted, provided below again, the '--set-filtered-gt-to-nocall' argument is supplied before '--max-nocall-number 42' so I was imagining the arguments will be dealt in the order they are listed.

    gatk SelectVariants \
    -R ~/reference_DM_PM4.03_G3_all_bowtie2/DM_v4.03_G3_all.fasta \
    -V $1.AllChr.final.fltd.variants.QD10.SOR5.MQ20.DP20t1000.GQ20.vcf.gz \
    --exclude-filtered \
    --select-type-to-include SNP \
    --restrict-alleles-to BIALLELIC \
    --set-filtered-gt-to-nocall \
    --max-nocall-number 42 \
    --selectExpressions 'vc.getGenotype("GPP116880").isHomRef()' \
    --exclude-non-variants \
    --add-output-vcf-command-line \
    -O $1.AllChr.final.fltd.slctd.variants.QD10.SOR5.MQ20.DP20t1000.GQ20.noCallGeno42.vcf.gz

    (2) I will try the fix you have suggested but I have some concerns with the no '2' suggestion you mentioned above. So if I use use '--max-nocall-number 21 --max-filtered-genotypes 21' and then any SNP with more than '21' for either categories will be eliminated even if the count for the other category is within acceptable limits. For example, I hope I am getting it right, if there are 25 nocall samples and 0 filtered-genotypes samples then the SNP in question will be excluded despite having far less nocall samples than the set figure of '42'

    (3) It looks the only option will be to process the vcf in two steps through the same SelectVariants tool, first implement '--set-filtered-gt-to-nocall' option and then process the output vcf through SelectVariants again with the remaining options OR apply all but '--max-nocall-number 42' option in the first step and the latter in the second step.

    As long as the issue just relates to the options' processing order then it is perfectly fine to work around this issue by any of the above approaches, I was more worried if there was/were more serious issues affecting this tool.

    Regards,
    Sanjeev

  • sanjeevkshsanjeevksh Member

    Hi @bhanuGandham

    Thanks for your help in resolving this matter,

    Regards,
    Sanjeev

Sign In or Register to comment.