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.

Why are some indels being called when in the read mapping there are none?

rafarios50rafarios50 ColombiaMember
edited May 2017 in Ask the GATK team

Greetings,

I have found some INDELS called that are not present in the read mapping and dont understand why it might happen or how to solve it.

I have tried to the calling with and without indel recalibration of the reads, but in both cases the same behaviour occurs.

The commands I have used are:
Without recalibration:

java -jar /source_installers/GenomeAnalysisTK.jar -T HaplotypeCaller -R ../../SNPS/503.fasta -I EfmD8t72_hdr.bam -o EfmD8t72_gatk_norecal.vcf

With recalibration:

java -Xmx2g -jar /source_installers/GenomeAnalysisTK.jar -T RealignerTargetCreator -R ../../SNPS/503.fasta -I EfmD8t72_hdr.bam -o EfmD8t72.intervals
java -Xmx4g -jar /source_installers/GenomeAnalysisTK.jar -T IndelRealigner -R ../../SNPS/503.fasta -I EfmD8t72_hdr.bam -targetIntervals EfmD8t72.intervals -o EfmD8t72_realigned.bam
java -jar /source_installers/GenomeAnalysisTK.jar -T HaplotypeCaller -R ../../SNPS/503.fasta -I EfmD8t72_realigned.bam -o EfmD8t72_gatk.vcf

The results I mentioned are:

gi|403040827|gb|AMBN01000124.1| 6975 . C CAG 7789.73 . AC=2;AF=1.00;AN=2;BaseQRankSum=-0.598;ClippingRankSum=0.000;DP=197;ExcessHet=3.0103;FS=0.000;MLEAC=2;MLEAF=1.00;MQ=58.04;MQRankSum=0.000;QD=32.80;ReadPosRankSum=1.735;SOR=0.313 GT:AD:DP:GQ:PL 1/1:1,175:176:99:7827,485,0
gi|403040827|gb|AMBN01000124.1| 6977 . ACC A 7780.73 . AC=2;AF=1.00;AN=2;BaseQRankSum=-0.567;ClippingRankSum=0.000;DP=197;ExcessHet=3.0103;FS=0.000;MLEAC=2;MLEAF=1.00;MQ=58.04;MQRankSum=0.000;QD=28.26;ReadPosRankSum=1.706;SOR=0.313 GT:AD:DP:GQ:PL 1/1:1,175:176:99:7818,485,0

The mapping of the reads as shown by CLC Genomics:
image

Also similar results happened in other regions of the same mapping.

Best Answer

  • SheilaSheila Broad Institute admin
    Accepted Answer

    @rafarios50
    Hi,

    Yes, you can do that. However, the bamout generation is quite time consuming (we only recommend using it for debugging purposes). You may consider only using it in areas where the regions are really messy.

    -Sheila

Answers

  • SheilaSheila Broad InstituteMember, Broadie admin

    @rafarios50
    Hi,

    Have you looked at the bamout?

    Thanks,
    Sheila

  • rafarios50rafarios50 ColombiaMember

    I'm not sure if I understand well, but the bamout would have the mapping of the reads right?
    Which would be the one obtained without realignment directly from BWA, in this case, or the obtained after the realignment, which is the one I'm showing in the image.

  • rafarios50rafarios50 ColombiaMember

    Sorry about the previous comment, you are right about the bamout and I could see the insertion and the deletion called from that alignment.

    If I wanted to obtain a consensus of variants called using other software and GATK this bamout should be the one used by the other callers aswell right?

  • SheilaSheila Broad InstituteMember, Broadie admin

    @rafarios50
    Hi,

    I am not sure what you mean? If you want to visualize the variants, then yes, you should use the bamout for that site because that is the alignment HaplotypeCaller uses to call variants.

    -Sheila

  • rafarios50rafarios50 ColombiaMember

    Well, to have a higher confidence on the variants, we are trying to obtain a consensus of variants called by different programs.
    Each one of those callers require a read mapping to do the calling, right?
    So my question would be: should I use that bamout file as input for the other variant callers rather than the bam file from BWA?

  • SheilaSheila Broad InstituteMember, Broadie admin
    Accepted Answer

    @rafarios50
    Hi,

    Yes, you can do that. However, the bamout generation is quite time consuming (we only recommend using it for debugging purposes). You may consider only using it in areas where the regions are really messy.

    -Sheila

Sign In or Register to comment.