A Sanger verified variant is not called by HaplotypCaller

Hello,

We have a trio which all three members were whole exome sequenced. We processed the samples by the standard pipeline (bwa mem alignment, picard mark duplicates and gatk joint variant calling by HaplotypeCaller). The version of GATK we used is v3.8.1. We can see there is a variant both from mother and the proband from IGV (see here) and the variant is Sanger verified both from proband and the mother. The position is well covered in both proband bam(94x from vcf FORMAT) and mother bam (207x from vcf FORMAT). However, in the variant file, the variant is called as de novo, so there is no variant called from mother.

17 71189251 . C G 1581.68 . AC=1;AF=0.167;AN=6;BaseQRankSum=1.59;ClippingRankSum=0.00;DP=343;ExcessHet=3.0103;FS=2.970;MLEAC=1;MLEAF=0.167;MQ=60.00;MQRankSum=0.00;QD=16.83;ReadPosRankSum=1.27;SOR=1.095 GT:AD:DP:GQ:PGT:PID:PL 0/1:46,48:94:99:0|1:71189182_A_G:1610,0,1851 0/0:41,0:41:99:.:.:0,99,1485 0/0:207,0:207:0:.:.:0,0,117

I have both bam file and bamoutput file from HaplotypeCaller sliced on this variant. I could upload if you tell me where to upload.

Could you help us to check what is the reason for the caller to miss this variant? Thanks.

Comments

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @ying_sheng_1
    Hi,

    I am assuming the last sample (3rd) is the mother? It looks like the PLs are 0,0 for the hom-ref and het genotypes, meaning the tool was not sure if the site is hom-ref or het (both are equally likely). You can try running the Genotype Refinement workflow to see if this changes. Also, can you check the base qualities at the site? Are they all pretty high? Can you also post the bamout with ~300 bases before and after the site of interest?

    Thanks,
    Sheila

  • ying_sheng_1ying_sheng_1 Member ✭✭
    edited May 2018

    Hi Sheila,

    Yes, you are right, the last sample is the mother. I attached the bamout with 300 bp up and downstream of the variant for all three members. So you have a positive case - proband, a negative case - father (also Sanger verified) and a problematic case - mother. I am trying the Genotype Refinement workflow now and will post the results later. For the base qualities, I think it is OK. I viewed on IGV and it looks fine.

    For the uploaded files, it seems the forum doesn't 'like' bam file, so I changed the file extension to "txt". When you use the files, you could remove ".txt" for each file.

    Thanks,

    Ying

  • ying_sheng_1ying_sheng_1 Member ✭✭

    The command to generate the bamout file is:
    java -Xmx2g -jar GenomeAnalysisTK.jar
    -T HaplotypeCaller
    -R human_g1k_v37_decoy.fasta
    -I bqsr.bam
    --emitRefConfidence GVCF
    --dbsnp dbsnp_137.b37.vcf
    --bamOutput proband.71189251.haplotype.bam
    --bamWriterType CALLED_HAPLOTYPES
    -L 17:71188951-71189551
    -o temp.g.vcf
    -forceActive
    -disableOptimizations

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @ying_sheng_1
    Hi Ying,

    Can you post IGV screenshots? I don't need the BAM files quite yet, I just need to see what the regions looks like in the BAMs and bamouts.

    Thanks,
    Sheila

  • hugueshugues NorwayMember

    My colleague Ying is away, but I could get you a screenshot:

  • ying_sheng_1ying_sheng_1 Member ✭✭
    edited May 2018

    My colleage is uploaded an IGV image from normal bam files for the variant. Now I loaded the IGV image from bamout output bams for the variant (the variant is in the middle).alt text

  • ying_sheng_1ying_sheng_1 Member ✭✭

    Hi Sheila,

    I also tried the genotype refinement you recommended. I tried "Derive posterior probabilities of genotypes" (CalculateGenotypePosteriors) on both joint variant calling vcf and single sample vcf just for mother.

    The command for joint variant calling vcf file is:

    java -Xmx2g -jar GenomeAnalysisTK.jar 
    -T CalculateGenotypePosteriors
    -R human_g1k_v37_decoy.fasta
    -V jointVariantCalling.vcf
    -ped jointVariantCalling.ped
    -o jointVariantCalling.withPosteriors.vcf
    

    The genotype for mother is changed to heterozygosity.

    The command for single sample vcf file is:

    java -Xmx2g -jar GenomeAnalysisTK.jar 
    -T CalculateGenotypePosteriors
    -R human_g1k_v37_decoy.fasta
    -V mother.vcf
    -o mother.withPosteriors.vcf
    

    The variant is still not called. I think it is because the variant is not in the input vcf file.

    Although the variant is found after CalculateGenotypePosteriors, we still wants to know the reason that this variant was not called from the beginning. Because we use the same way to do variant calling on single sample, and it seems that even with CalculateGenotypePosteriors, the variant in single sample is not recovered.

  • ying_sheng_1ying_sheng_1 Member ✭✭

    Hi Sheila,

    Sorry for pushing this, is there any progress on this case?

    thanks!

    Ying

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin
    edited May 2018

    @ying_sheng_1
    Hi Ying,

    I don't think the Genotype Refinement workflow will work well on a single sample unless you input a population resource VCF. It worked on the trio because you input a ped file which helps inform the priors and the re-calculation of the PLs.

    In this case, when you are really interested in correct genotypes, and not just variant sites, it would be best to run the Genotype Refinement workflow. Unfortunately, in some cases like the mother, there is just not enough clear evidence for the tool to make a proper genotype call.

    -Sheila

  • ying_sheng_1ying_sheng_1 Member ✭✭

    Hi Sheila,

    For this case, I have some updates. I have rerun variant calling with "samtools" and "UnifiedGenotyper" on the variant site, the mother actually got the right genotype (0/1), and the GQ from UnifiedGenotyper is 99.

    I have rerun the Genotype Refinement workflow on the mother sample with a supporting file from 1000G (ALL.wgs.phase3_shapeit2_mvncall_integrated_v5b.20130502.sites.vcf.gz). The variant is not recovered. I don't think it will be recovered if the input vcf doesn't contain the variant.

    I also got the pileup on the position from the mother bamout file and check the base quality. I didn't see the differences on the base qualities between those contains the reference allele and those with alternative allele, both have base qualities over than 20. The mapping qualities from the bam (reads in +/- 300bp of the region) I sent to you is also good.

    I could not find a reason to explain why HaplotypeCaller could not decide the genotype.

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @ying_sheng_1
    Hi Ying,

    I see. Thanks for doing some extra digging. It sounds like this is a reassembly issue. My colleague is interested in improving the reassembly algorithm, so I will take a look at these BAMs you provided and see if they can help him.

    -Sheila

  • ying_sheng_1ying_sheng_1 Member ✭✭

    Thanks, Sheila!

  • ying_sheng_1ying_sheng_1 Member ✭✭

    Hi Sheila,

    Is there any progress for this case?

    thanks,

    Ying

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    Hi @ying_sheng_1, sorry for the delayed response. There is nothing we can do immediately to resolve this case; however the developer who is responsible for the assembly code implicated here has been collecting multiple cases like this where the program performs poorly, and he plans to use them to improve that functionality. Unfortunately I can't give you a timeframe for when the improvements might become available as it's difficult to predict how long it will take to identify the source of the problem, then devise a solution, implement and test it thoroughly. It can take anywhere between a few weeks and several months.

  • ying_sheng_1ying_sheng_1 Member ✭✭

    Hi Geraldine,

    Thanks for explaining this. Hope this will help to improve the method.

    nice summer,

    Ying

Sign In or Register to comment.