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.

Reference calls absent at certain locus for half of the samples with >500 coverage

Hi @bhanuGandham

I have now uploaded a new snippet file (snippet_jayaramv_Aug13) with a corrected 80019_snippet.g.vcf, merged gvcf and joint called vcf record.

1) 80019_snippet.g.vcf is an empty file in the bug report.

Reuploaded correct version

2) I looked at the bamout for the other 6 samples, and all samples have no reads covering 14 95590730 region. That is the reason for the no call.

I know, I didn't upload that region for the 7 samples. The region is 16:50824000-50829000 .

3) And for this region:16:50826538, can you please post the joint called variant vcf record.

Done


Just to reiterate my previous question in a better way

Study design
GATK version 4 following best practice.
600 samples of targeted gene panel using nextera (100 genes)
I used join calling after merging gvcfs.

Looking across the 600 samples [for example at variant 16:50826538, I have many similar loci across the panel] I notice the following:

The final vcf for some variants around 30% samples have no calls (./.) and the rest are (0/0) or (1/0) if there are any...

The regions with (./.) do not appear to have poor read depths and looks quite similar in quality to the samples that are called (0/0).

So my question is why do a large percentage of samples have no calls at certain variant locus for no apparent reason?

I have uploaded a 5 Kb region for 3 samples with no calls, 3 samples with reference calls and one sample with a heterozygote call with bam, gvcf , bamout, merged gvcf and vcf files for each.

-L 16:50824000-50829000 was taken as an example.
This region is just one of the many regions of my targeted resequencing just to highlight my issue.


I have uploaded a 5kb region for seven samples which show this problem clearly.


All samples have good read depths but three samples have no calls at 16:50826538, but three are 0/0 and one sample is (0/1).

16:50826538
79880: ./.
80019: ./.
80020: ./.

80008: 0/0
80009: 0/0
80010: 0/0

80600: 0/1

My question why do some samples have no calls while others have.

The VCF is before the step "Hard-filtering germline short variants" which may remove this variant because of QD<2. However I have many other variants that have say for example QD>20 that have the same no calls issue.


Thanks again for looking into this issue, I really appreciate it and I am struggling to understand why this is happening at various different loci

Answers

  • jayaramvjayaramv Member
    Hi @bhanuGandham

    Thanks for your prompt reply.

    I did go through the bamout files and visualise them in the IGV viewer before I uploaded the data to ftp site. I have already noticed that there are not active regions in that site for some samples :( .

    1. I appreciate that 79880,80019 and 80020 dont have any reads in the bamout that align to that position but so does 80008,80009 and 80010. However those three have good Phred quality scores and call the reference allele as 0/0. Why is this happening?

    2. If you look at the bam files themselves of all 7 samples you can see that the region has great coverage in all samples and the quality looks identical.


    So in spite of no reads in the bamout stage samples still give reference calls. So can you kindly have a look at the bam files themselves and comment why there are no reference calls at this site.


    Apologies for being a pain.
  • bhanuGandhambhanuGandham Cambridge MAMember, Administrator, Broadie, Moderator admin

    Hi @jayaramv

    I looked at your data. If you look at the bamout for the samples 79880, 80019 and 80020 at the 16:50826538 region you will see that no reads align to this position. This is the reason for the no calls. And if you looking for why you see no reads in the bamout, that description is given in these docs:
    https://software.broadinstitute.org/gatk/documentation/article?id=11068
    https://gatkforums.broadinstitute.org/gatk/discussion/11076/local-re-assembly-and-haplotype-determination-haplotypecaller-mutect2#latest

  • bhanuGandhambhanuGandham Cambridge MAMember, Administrator, Broadie, Moderator admin

    Hi @jayaramv

    The docs I shared will explain why and how haplotypecaller realigns reads in the active regions. Take a look at the docs that should explain it.

  • bhanuGandhambhanuGandham Cambridge MAMember, Administrator, Broadie, Moderator admin

    Also, can you please post the exact command you are using to create the gvcf files.

  • jayaramvjayaramv Member
    HI @bhanuGandham

    thanks for the reply and the links. I understand why there are no active regions for the reference calls and that's why there is nothing in the bamout. The question I have is why these samples are ./.

    79880: ./.
    80019: ./.
    80020: ./.

    and these are 0/0

    80008: 0/0
    80009: 0/0
    80010: 0/0

    Looking at the bam files they appear to both have high depth. The full call in the vcf for for 79880 which is ./. is : GT:AD:DP:GQ:PL ./.:625,0:625:.:0,0,0

    and for 80008 which is 0/0 is 0/0:151,0:151:99:0,112,1800

    I know HaplotypeCaller assess a call/no call on many factors but it's not clear why it doesn't call a reference in these cases. There are several examples of calls like this in the vcf. The reason why I am worried is that we will be conducting a case control analysis and the call rate will affect the variant frequency and our confidence in the heterozygous calls. Sufficient no calls will also mean that we reject a site (e.g. < 85% callrate).

    Here are the commands (all GATK4.1)

    gatk --java-options -Xmx8G HaplotypeCaller -R GATK_b37_reference/b37/human_g1k_v37_decoy.fasta -I 80019_snippet.bam -D GATK_b37_reference/b37/dbsnp_138.b37.vcf -O 80019_snippet.g.vcf -ERC GVCF -bamout 80019_snippet_bamout.bam -L 16:50824000-50829000

    gatk GenotypeGVCFs -R GATK_b37_reference/b37/human_g1k_v37_decoy.fasta -D GATK_b37_reference/b37/dbsnp_138.b37.vcf --variant snippiet.g.vcf -L 16:50824000-50829000

    gatk CombineGVCFs -R GATK_b37_reference/b37/human_g1k_v37_decoy.fasta -O snippiet.g.vcf --variant 79880_snippet.g.vcf --variant 80020_snippet.g.vcf --variant 80008_snippet.g.vcf --variant 80009_snippet.g.vcf --variant 80010_snippet.g.vcf --variant 80600_snippet.g.vcf --variant 80019_snippet.g.vcf -O snippet.vcf

    best,
    Jay
  • jayaramvjayaramv Member
    Hi @bhanuGandham ,
    Just wondering if you had any insights to the problem.
    Best, Jay
  • bhanuGandhambhanuGandham Cambridge MAMember, Administrator, Broadie, Moderator admin

    Hi @jayaramv

    I looked at your data. If you look at the bamout for the samples 79880, 80019 and 80020 at the 16:50826538 region you will see that no reads align to this position. This is the reason for the no calls. And if you looking for why you see no reads in the bamout, that description is given in these docs:
    https://software.broadinstitute.org/gatk/documentation/article?id=11068
    https://gatkforums.broadinstitute.org/gatk/discussion/11076/local-re-assembly-and-haplotype-determination-haplotypecaller-mutect2#latest

Sign In or Register to comment.