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.
Attention:
We will be out of the office for a Broad Institute event from Dec 10th to Dec 11th 2019. We will be back to monitor the GATK forum on Dec 12th 2019. In the meantime we encourage you to help out other community members with their queries.
Thank you for your patience!

Variant calling: high (and strange) number of alternative allele

Deat GATK team,

I am calling variant on a trio (mother, father and offspring) of Macaca mulatta. I have whole genome sequencing 60X for each individual. I use GATK 4.0.7.0, I call variant with HaplotypeCaller BP-RESOLUTION, combine with GenomicDBimport per chromosomes and genotype with GenotypeGVCF.

I am interested in the number of sites where I have only reference allele (AD=0 for the alternative) and the number of sites where I have some reads supporting ALT allele (AD > 0) in the parents.

I found a lot of sites (for each individuals) where I have AD>0 in the gvcf file (per indiviuals, the combined one and after genotyping). I looked at each site that are HomRef and for each individuals less than 30% of the HomRef sites have AD=0 for the alternative allele. I know that HaplotypeCaller does a realignement step that may change the positions of the reads, but 70% of the sites that have AD>0 seems a lot. I looked back at the BAM file and those alternative alleles don’t seem to be there. I try to call again using the bam.out option, and here again I don’t see so many alternative alleles. However, I see that sometimes on a read where there were no alternative allele on the bam input there is an alternative allele on the output.
Also I have tried samtools mpileup and in this case almost 90% of the HomRef sites are AD=0 for alternative allele.

Just as an example bellow is the VCF output from HaplotypeCaller for one individual and then there is a picture of both the input bam file and the output bam file.
For chr1 pos 24203380 the ref is A and I have:
Vcf --> DP=96 AD=92,4
Bam input --> DP 93, 92,1 (N)
Bam out --> DP=80, 79,1 (N)

chr1 24203380 . A <NON_REF> . . . GT:AD:DP:GQ:PL 0/0:92,4:96:57:0,57,5771 chr1 24203381 . G <NON_REF> . . . GT:AD:DP:GQ:PL 0/0:90,5:95:0:0,0,5897 chr1 24203382 . C <NON_REF> . . . GT:AD:DP:GQ:PL 0/0:92,3:95:78:0,78,6075 chr1 24203383 . A <NON_REF> . . . GT:AD:DP:GQ:PL 0/0:92,3:95:68:0,68,6127

Just in case here is my code:
gatk --java-options "-XX:ParallelGCThreads=16 -Xmx64g" HaplotypeCaller -R /PATH/rheMac8.fa -I /PATH/R01068_sorted.merged.addg.uniq.rmdup.bam -O /PATH/R01068_res.g.vcf -ERC BP_RESOLUTION \

I don’t know why I have this high number of alternative alleles and how to get read of them to have the 'real' number of alternative allele per position. The problem persists on the genotyping vcf files with some alternative alleles that are not present on any bam (input or HaplotypeCaller output).

I hope I gave you enough details so you have a clear idea of my problem and will be able to help me.
Best,

Best Answer

Answers

  • AdelaideRAdelaideR Member admin

    Hello @Lucie_B - The HaplotypeCaller is a good first pass method to use to find your variants. You might want to consider adding the Genotype Refinement Workflow to this process to make sure that your child genotype reflects the parental probabilities of observing an allele.

    For a sample for which both parent genotypes are available, the child’s genotype can be supported or invalidated by the parents’ genotypes based on Mendel’s laws of allele transmission. Even the confidence of the parents’ genotypes can be recalibrated, such as in cases where the genotypes output by HaplotypeCaller are apparent Mendelian violations.

    This should probably remove some of the alternative alleles currently seen in your analysis.

    Another procedure to consider is adjusting your hard filtering options.

    This will reduce the number of calls to those that fall within a certain threshold value.

  • Lucie_BLucie_B Member

    Dear @AdelaideR

    Thanks for the advice of using the Genotype Refinement Workflow. I have looked at it and it does reduce the number of Mendelian Violation seen in my dataset.

    However I am not sure this will be helpful on my problem. My issue is mostly on the sites for which all samples are HomRef. I have a non normal proportion of those sites with at least one alternative alleles seen. See below an example:

    chr1 34 . T <NON_REF> . . . GT:AD:DP:GQ:PL ./.:88,1:89:99:0,120,1800 ./.:65,2:67:65:0,65,4228 ./.:72,1:73:99:0,120,1800
    chr1 35 . T <NON_REF> . . . GT:AD:DP:GQ:PL ./.:86,1:87:99:0,120,1800 ./.:63,3:66:0:0,0,3997 ./.:72,1:73:99:0,120,1800
    chr1 36 . T <NON_REF> . . . GT:AD:DP:GQ:PL ./.:84,2:86:99:0,113,5502 ./.:63,2:65:58:0,58,4124 ./.:74,3:77:34:0,34,4380
    chr1 37 . A <NON_REF> . . . GT:AD:DP:GQ:PL ./.:87,1:88:99:0,120,1800 ./.:62,2:64:60:0,60,4034 ./.:77,1:78:99:0,120,1800
    chr1 38 . G <NON_REF> . . . GT:AD:DP:GQ:PL ./.:87,1:88:99:0,120,1800 ./.:63,3:66:0:0,0,3920 ./.:79,1:80:99:0,120,1800
    chr1 39 . G <NON_REF> . . . GT:AD:DP:GQ:PL ./.:90,1:91:99:0,120,1800 ./.:61,3:64:0:0,0,3895 ./.:77,2:79:91:0,91,4600
    chr1 40 . A <NON_REF> . . . GT:AD:DP:GQ:PL ./.:87,2:89:99:0,120,1800 ./.:57,4:61:0:0,0,3493 ./.:77,3:80:53:0,53,4302

    All those sites are HomRef but with few alternative alleles that I don't see on the Bam file (neither the input ones or the output ones). I have more than 70% of the sites with at least one alternative allele, but only 10 % of the sites have alternative allele in the Bam files.

    Would you have any idea where those alternatives allele could come from? As they are neither in the bam in or bam out I find it very strange. I wonder if this is normal or if something could have went wrong during variant calling.
    Best,

  • AdelaideRAdelaideR Member admin

    @Lucie_B

    This is interesting behavior. I am going to do a little more research.

    My first hunch is that this is happening during the realignment phase of the HaplotypeCaller.

    It would be helpful to see a screenshot of the bams to make sure that the alternative alleles are not present.

    I will keep doing some research in the meantime.

  • Lucie_BLucie_B Member

    Dear GATK team,

    I come back to you has I can not fix my problem and feel very stuck here.
    I still have this issue with a lot of sites (more than 80% of the whole genome) with at least one alternative allele. Link to this I saw that many sites have a GQ=0, while those sites seem to be HomRef and with a very high coverage. This problem occurs always in blocks of 3 to 7 sites following each other, but when I look at the bam file I don't see what could have causes this GQ of 0.

    Here is an exemple, this is the gvcf file in bp resolution mode for one sample:

    Position 150 has a good GQ while the following positions are GQ=0

    Here is what the bam file looks like in IGV:

    All the reads seem too have a good mapping quality.

    Would you have any idea of why this is happening?
    Best regards,

    Lucie

  • bshifawbshifaw Member, Broadie, Moderator admin

    @Lucie_B

    According to the posts above we'll need a snippet of the bam and commands used to better help you. The instructions in this article should get you started in sharing the necessary material.

Sign In or Register to comment.