Heads up:
We’re moving the GATK website, docs and forum to a new platform. Read the full story and breakdown of key changes on this blog.
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.

BAM with hight depth : variant called as HOM_VAR but there are almost no ALT read

lindenblindenb FranceMember ✭✭

Hi GATK Team,

I'm using 4.1.2.0 / HaplotypeCaller with the following command with a Haloplex BAM (high depth)

 gatk --java-options "-Xmx5g -Djava.io.tmpdir=." HaplotypeCaller --minimum-mapping-quality "20" -L "15:42170528-42172528" -R human_g1k_v37.fasta -I jeter.bam -O jeter.vcf  --dont-use-soft-clipped-bases  --max-reads-per-alignment-start 1000

the generated vcf contains a HOM_VAR variant at 15:42171528

15  42171528    .   G   T   2573.03 .   AC=2;AF=1.00;AN=2;BaseQRankSum=-5.186;DP=971;ExcessHet=3.0103;FS=9.978;MLEAC=2;MLEAF=1.00;MQ=60.00;MQRankSum=0.000;QD=3.70;ReadPosRankSum=-7.353;SOR=7.656  GT:AD:DP:GQ:PL  1/1:71,624:695:99:2587,626,0

while I can mostly see REF alleles:

$  samtools mpileup -f human_g1k_v37.fasta jeter.bam | grep 42171528

15  42171528    G   1060    ,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,t,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,t,,t,,c,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,................................................................,    B>?B???>???:>?>?BB:8>><???:???>??B????????????????????:?B???????:?????????>[email protected][email protected]:?>@[email protected]>??>:?????????>>[email protected]???>????<B????:>@[email protected]??:???????<[email protected]???<????????????>@:??1B?<A?>[email protected][email protected][email protected]???6?1???>??B??????????????????????:??????B:??>??:???????8?:[email protected]??5B5???????8?:>?????><8?B8?:[email protected]??>??<[email protected][email protected]?<[email protected]??B?????><[email protected]>?>??>??>[email protected]???/?:[email protected][email protected]?>??>?:8B8??????>?:?B????????BBB??????????>[email protected]>???:?????>@?????????????B?:?<>?????>?>???:[email protected]???<?????????:?>3??>??:?>>?8>???:???????????:???:[email protected]???????????>???????????5????????8???B??????????>[email protected][email protected]>;A?88A;@[email protected]>=7BA?AA>[email protected][email protected]@[email protected]@??;[email protected];[email protected]?BAAA?ABAA?AAA?AA???D?>A?AAADA8BAAAD?AA?AAD?>AAD>A/[email protected];AABA8AAAAA7A?A?AA???;[email protected]>AA;[email protected]@;[email protected][email protected]?A?AAA>[email protected]@[email protected]@AAA?AA6ABAB>[email protected]?ADAA??A??AA2;[email protected][email protected][email protected]?A???A?A;;[email protected][email protected]=AAAAAAAAAA??>@[email protected]/?A;[email protected]@@[email protected][email protected]?DAAA>[email protected];;9ADD/:B>6>8;>:6:@8;/[email protected]@[email protected]^cgggh`ge__cgg^eegiggZegkgeggihgaggggjgVgegggeiXgeigX_ggegdajB

Am I doing anything wrong ? If it helps, I put some minimal files here: https://nextcloud-bird.univ-nantes.fr/index.php/s/HYsx26oRmrJNwAd

Best Answers

  • lindenblindenb France ✭✭
    Accepted Answer

    @bshifaw
    thank you very much for your answer.

    1) and 2) as stated above, it's haloplex , it always looks like this because it's a specific capture using enzymes: it always looks like a set of cluster of reads.
    3) this is a subset of our BAM, I've narrowed the bam to a region around the problematic SNP.

    I'm not sure I understand all your description of the graph , using --adaptive-pruning seems to work: the HOM_VAR variant was removed.

    Thank you so much !
    P.

Answers

  • bshifawbshifaw Member, Broadie, Moderator admin

    @lindenb
    Thanks for the files, would you also please provide a bamout using the following instructions Article5484

  • lindenblindenb FranceMember ✭✭

    @bshifaw thanks for your help. I added the file jeter.bamout.bam and its' index in my shared directory.

  • lindenblindenb FranceMember ✭✭

    @bshifaw ops sorry, i didn't use -forceActive -disableOptimizations . give me 5 minutes...

  • lindenblindenb FranceMember ✭✭

    @bshifaw ok done with --force-active --disable-optimizations

  • lindenblindenb FranceMember ✭✭

    @bshifaw just a friendly ping. Do you have any news for this problem ? Thank you.

  • bhanuGandhambhanuGandham Cambridge MAMember, Administrator, Broadie, Moderator admin

    Hi @lindenb

    @bshifaw is looking into this issue and will get back to you shortly.

  • bshifawbshifaw Member, Broadie, Moderator admin

    Haplotypecaller performs a local assembly within active regions to improve calling, so its best to look at the bamout file instead of the input file in IGV. The bamout file is the intermediate file after HC performs local assembly and is what HC uses when calling variants. Looking at the bamout in IGV it looks like all of the reads are reference and haplotypecaller is still calling a variant at this location. I've asked the dev team and they mentioned if the site is

    perfectly phased then it might be possible under some circumstances.

  • lindenblindenb FranceMember ✭✭

    @bshifaw thank you for your answer but I'm afraid it doesn't explain why I get a variant.

    As far as I understand, in your screenshot you're showing the stats from the VCF file which is HOM_VAR (my problem).

    But if you're looking at the bamout, after the local assembly, in IGV, you'll see that** 98% of the reads carry the REF allele (G).**

    Why this HOM_VAR mutation ?

  • lindenblindenb FranceMember ✭✭

    @bshifaw even with samtools mpilup, i can see all those REF='G' in the bamout file:

    $ samtools mpileup jeter.bamout.bam  | grep 42171528
    [mpileup] 2 samples in 1 input files
    15  42171528    N   972 ttttggttggggggtttgtgtggtgggtggggtgggggtgtttgggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggtgtgggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggggtgggggggggggggggggGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGggggggggggggggggggggggggg
    
  • bshifawbshifaw Member, Broadie, Moderator admin

    I'll see if I can dig up some more info from the dev team.

  • bshifawbshifaw Member, Broadie, Moderator admin

    @lindenb ,

    We're still looking into your results, would you mind rerunning your command with the following parameters added
    --graph-output
    --debug-graph-transformations

    --graph-output should point to a file with the .dot extension
    --debug-graph-transformations - should do this for a subset just corresponding to the active region in question, then it should print the graph to where the user specifies.

    Post the .dot file produced

  • lindenblindenb FranceMember ✭✭

    thanks; I added --graph-output jeter.dot --debug-graph-transformations to my command. I added all the generated *.dot files to my shared directory.

    thanks again for your help.

  • lindenblindenb FranceMember ✭✭
    Accepted Answer

    @bshifaw
    thank you very much for your answer.

    1) and 2) as stated above, it's haloplex , it always looks like this because it's a specific capture using enzymes: it always looks like a set of cluster of reads.
    3) this is a subset of our BAM, I've narrowed the bam to a region around the problematic SNP.

    I'm not sure I understand all your description of the graph , using --adaptive-pruning seems to work: the HOM_VAR variant was removed.

    Thank you so much !
    P.

Sign In or Register to comment.