Notice:
If you happen to see a question you know the answer to, please do chime in and help your fellow community members. We encourage our fourm members to be more involved, jump in and help out your fellow researchers with their questions. GATK forum is a community forum and helping each other with using GATK tools and research is the cornerstone of our success as a genomics research community.We appreciate your help!

Test-drive the GATK tools and Best Practices pipelines on Terra


Check out this blog post to learn how you can get started with GATK and try out the pipelines in preconfigured workspaces (with a user-friendly interface!) without having to install anything.
Attention:
We will be out of the office on November 11th and 13th 2019, due to the U.S. holiday(Veteran's day) and due to a team event(Nov 13th). We will return to monitoring the GATK forum on November 12th and 14th respectively. Thank you for your patience.

HaplotypeCaller not calling SNP at position that appears should be called.

Using the following line I'm calling SNPs.
java -Xmx4g -jar ${GATKPath} -R $ref -T HaplotypeCaller -ploidy 1 -I downsample.bam -o hapreadyAll.vcf -bamout bamout.bam -dontUseSoftClippedBases

Shown in the screen shot it seems for this regions a bamout file along with a SNP should be output. Is here a HC command to force this region and others like it to output a SNP call?

Answers

  • pdexheimerpdexheimer Member ✭✭✭✭

    If you're not getting a bamout at that position, that means that HC believes strongly enough that the region is invariant that it skips processing it altogether. There are several levels of optimization, and flags to disable each of them. If you run with --disableOptimizations --dontTrimActiveAlignments --forceActive you'll disable all of them and should get something in the bamout. These will drastically slow down execution, you should probably run it only over the ~1kb or so surrounding your site (don't go too small, HC needs to see some context)

  • tstubertstuber Member

    You are right, very slow. Over 90 minutes to get 8kb ran.

    SNP still not called with the following additions. Any other suggestions?

    --disableOptimizations --dontTrimActiveRegions --forceActive

  • pdexheimerpdexheimer Member ✭✭✭✭

    It shouldn't change the calls, it should change the contents of the bamout file, which should help to diagnose why nothing was called

  • tstubertstuber Member

    Update: There was a bamout file but no SNP called.

  • SheilaSheila Broad InstituteMember, Broadie admin

    @tstuber
    Hi,

    Can you post the bam and bamout file screenshots with the extra arguments?
    Also, can you check that the base qualities and mapping qualities are good in the region?

    Thanks,
    Sheila

  • tstubertstuber Member

    Attached are the bam and bamout files. See SNP at position 9. Overall mapping qualities are good.

    Thanks,
    Tod

  • SheilaSheila Broad InstituteMember, Broadie admin

    @tstuber
    Hi,

    Is this the same region you posted in the beginning of the thread? What happens if you don't add -dontUseSoftClippedBases? It looks like all the SNP Ts are in the very beginning of reads which is not good.

    -Sheila

  • tstubertstuber Member

    Not the same region. Same issue though. I hope this is okay.

    I tried using just...
    java -Xmx4g -jar ${gatk} -R $ref -T HaplotypeCaller -I mybamfile.bam -o $n.hap.vcf -bamout $n.bamout.bam
    with all flags removed, but SNP was not called.

  • SheilaSheila Broad InstituteMember, Broadie admin

    @tstuber
    Hi,

    Yea, unfortunately, I don't think there is anything we can do to help you out here. You can try lowering the confidence thresholds, but it looks like Haplotype Caller just doesn't have enough evidence to call the SNP. Just as a check, when you run Haplotype Caller in BP_RESOLUTION mode or in GVCF mode, what exactly is shown at the sites? Are the PLs all 0s? If the PLs are all 0s, the default is to call a hom ref site.

    -Sheila

  • tstubertstuber Member

    Using haplotype gvcf
    segment3_PA_consensus-KF467227 9 . C <NON_REF> . . END=9 GT:DP:GQ:MIN_DP:PL 0/0:21:0:21:0,0,0

    Using haplotype bp_resolution
    segment3_PA_consensus-KF467227 9 . C <NON_REF> . . . GT:AD:DP:GQ:PL 0/0:0,21:21:0:0,0,0

    Using unifiedgenotyper emit all
    SeqPAcons 9 . C T 1646.77 . AC=2;AF=1.00;AN=2;DP=70;Dels=0.03;FS=0.000;HaplotypeScore=2.9277;MLEAC=2;MLEAF=1.00;MQ=60.00;MQ0=0;QD=23.53;SOR=1.263 GT:AD:DP:GQ:PL 1/1:0,68:68:99:1675,153,0

    I'm surprised the UG is calling it but the HC isn't.

    Glad to have help troubleshooting and trying all options.

    Thanks!

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    Out of curiosity, why is the contig name different between the HC and the UG runs? Are you running with a different reference file?

  • tstubertstuber Member

    The reference headers are different, but the sequences are the same.

    Thanks,

  • SheilaSheila Broad InstituteMember, Broadie admin

    @tstuber
    Hi,

    The fact that in BP_RESOLUTION mode, the AD shows 21 non-ref, but the PLs are all 0 suggests there is not good enough evidence to call the variant. Again, when PLs are all 0, the default is to call the site hom-ref. Maybe if you add a few more samples with a variant at the site, it will be called as a variant site.

    -Sheila

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    Be sure to check the base qualities at the site, not just the mapping qualities. HC imposes a stringent base quality filter, but the value can be lowered, iirc.

  • tstubertstuber Member

    Thank you, very good information.

Sign In or Register to comment.