Does the CalculateGenotypePosteriors change parents' genotype?

zwangzwang torontoMember

For example, after the gatk recalibration step, the genotype of a parent, which was 0/0, was changed to 0/1 after running CalculateGenotypePosteriors?

Answers

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @zwang
    Hi,

    Yes, it is possible for the genotype to change after CalculateGenotypePosteriors. Have a look at the documentation for more information.

    -Sheila

  • zwangzwang torontoMember

    We analyzed a trio following the Genotype Refinement workflow: https://www.broadinstitute.org/gatk/guide/article?id=4723.
    We found a weird record.

    Before applying CalculateGenotypePosteriors, the mother has genotype 0/0 (2nd to the last column)
    chr1 96925200 . A T 530.13 PASS AC=1;AF=0.167;AN=6;BaseQRankSum=0.599;ClippingRankSum=1.20;DP=87;FS=0.000;GQ_MEAN=486.00;MLEAC=1;MLEAF=0.167;MQ=60.00;MQ0=0;MQRankSum=-8.750e-01;NCC=0;QD=18.93;ReadPosRankSum=0.921;SOR=0.874;VQSLOD=16.04;culprit=MQ GT:AB:AD:DP:GQ:PL 0/1:0.460:13,15:28:99:561,0,486 0/0:.:25,0:25:34:0,34,969 0/0:.:31,0:31:68:0,68,1097

    After applying CalculateGenotypePosteriors, the mother's genotype is 0/1,
    chr1 96925200 . A T 530.13 PASS AC=2;AF=0.333;AN=6;BaseQRankSum=0.599;ClippingRankSum=1.20;DP=87;FS=0.000;GQ_MEAN=486.00;MLEAC=1;MLEAF=0.167;MQ=60.00;MQ0=0;MQRankSum=-8.750e-01;NCC=0;QD=18.93;ReadPosRankSum=0.921;SOR=0.874;VQSLOD=16.04;culprit=MQ GT:AB:AD:DP:GQ:JL:JP:PL:PP 0/1:0.460:13,15:28:99:0:24:561,0,486:527,0,545 0/1:.:25,0:25:25:0:24:0,34,969:25,0,935 0/0:.:31,0:31:34:0:24:0,68,1097:0,34,1063

    But we know, for sure, that the mother's genotype is homo reference, the kid has a de novo call. Here is the command I used for running the CalculateGenotypePosteriors

    java -Xmx5g -jar GenomeAnalysisTK.jar
    -T CalculateGenotypePosteriors
    -R ucsc_hg19.fa
    -V joint_recal.vcf
    --skipPopulationPriors
    -ped fam.ped
    -o joint_recal_gposteriors.vcf

    We used samtools mpileup to look at the reads at that location on both parents, there are no supporting reads for alternative allele.

    Any idea on what is going on?

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @zwang
    Hi,

    Can you tell me how you generated the trio VCF? Also, can you please post a screenshot of the bam file in that region? Please also post the bamout file in the region.

    Thanks,
    Sheila

  • zwangzwang torontoMember

    For each sample, here the command to call

    java -Xmx11G -Djava.io.tmpdir=/localhd/468805[9]/tmp.eWZ7LiizrR
    -jar GenomeAnalysisTK.jar
    -T HaplotypeCaller
    -R ucsc_hg19.fa
    -o child.chr8.htc.vcf.ungenotyped.VCF
    -GQB 10 -GQB 30 -stand_emit_conf 10 -A AlleleBalance -A AlleleBalanceBySample -A BaseCounts -A ClippingRankSumTest -A Coverage -A DepthPerAlleleBySample -A DepthPerSampleHC -A FisherStrand -A GCContent -A GenotypeSummaries -A HardyWeinberg -A HomopolymerRun -A InbreedingCoeff -A LowMQ -A MappingQualityRankSumTest -A MappingQualityZero -A MappingQualityZeroBySample -A NBaseCount -A QualByDepth -A RMSMappingQuality -A ReadPosRankSumTest -A SpanningDeletions -A StrandBiasBySample -A StrandOddsRatio -A TandemRepeatAnnotator -A VariantType
    --dbsnp dbsnp_138.ucsc.hg19.vcf --genotyping_mode DISCOVERY --emitRefConfidence GVCF --variant_index_type LINEAR --variant_index_parameter 128000 -L chr8 -nct 8 -I input.bam

    Since I made calls by chromosome, I next merged the calls.

    Then the joint genotyping is:
    java -Xmx5g -jar GenomeAnalysisTK.jar -R ucsc_hg19.fa -T GenotypeGVCFs -o joint.vcf --dbsnp dbsnp_138.ucsc.hg19.vcf -nt 8 mother.g.vcf --variant:VCF father.g.vcf --variant:VCF child.g.vcf

    I will generate the bamout file shortly.

  • zwangzwang torontoMember

    If we upload the bamout files, are they private? Basically, there was no reads from mother and father.

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    The uploads go to a private directory that can only be accessed by our staff. We can't 100% guarantee we'll never get hacked, but we do make every reasonable effort to prevent data leaks. And we do delete user data when we're done with it, unless we are explicitly given permission to use some of it for systematic testing purposes.

  • zwangzwang torontoMember

    Ok. How to upload the files?

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @zwang
    Hi,

    If you post the bam file screenshot to the forum, they are viewable by the public. Geraldine is referring the uploading to the FTP when we ask you to upload a bug report. For now, I am requesting a small screenshot of the bam file and bamout file of the region in question. I want to see why the call was changed from homozygous reference to heterozygous. To do that, I first need to see what the region looks like in the bam file. The GQ from HaplotypeCaller is 34 which is not very high, so I wanted to see what is going on in the region.

    How do you know for sure the mother is homozygous reference?

    Thanks,
    Sheila

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @zwang
    Hi again,

    To upload the files, you can click on "Attach a file" below the comment box. Then select the screenshots you would like to upload. We really prefer this way because it is easy for us to view, and the region is not easily identifiable. However, if you prefer to keep things private, you can submit the files to our FTP. Instructions are here. You don't have to submit all the files required in the instructions, just the relevant bam files.

    -Sheila

  • zwangzwang torontoMember

    Here is the screen shot, as I mentioned before there is no reads from mother. And one thing, when I used the -bamout option, I did not use -nct 8; the program says it does not support that option with bamout.

  • zwangzwang torontoMember

    Here is the problematic site: before applying CalculateGenotypePosteriors the mother's genotype (the bold) is 0/0, after CalculateGenotypePosteriors it became 0/1

    joint_recal.vcf:chr1 96925200 . A T 530.13 PASS AC=1;AF=0.167;AN=6;BaseQRankSum=0.599;ClippingRankSum=1.20;DP=87;FS=0.000;GQ_MEAN=486.00;MLEAC=1;MLEAF=0.167;MQ=60.00;MQ0=0;MQRankSum=-8.750e-01;NCC=0;QD=18.93;ReadPosRankSum=0.921;SOR=0.874;VQSLOD=15.95;culprit=MQ GT:AB:AD:DP:GQ:PL 0/1:0.460:13,15:28:99:561,0,486 0/0:.:25,0:25:34:0,34,969 0/0:.:31,0:31:68:0,68,1097

    joint_recal_gposteriors.vcf:chr1 96925200 . A T 530.13 PASS AC=2;AF=0.333;AN=6;BaseQRankSum=0.599;ClippingRankSum=1.20;DP=87;FS=0.000;GQ_MEAN=486.00;MLEAC=1;MLEAF=0.167;MQ=60.00;MQ0=0;MQRankSum=-8.750e-01;NCC=0;QD=18.93;ReadPosRankSum=0.921;SOR=0.874;VQSLOD=15.95;culprit=MQ GT:AB:AD:DP:GQ:JL:JP:PL:PP 0/1:0.460:13,15:28:99:0:24:561,0,486:527,0,545 0/1:.:25,0:25:25:0:24:0,34,969:25,0,935 0/0:.:31,0:31:34:0:24:0,68,1097:0,34,1063

    BTW, the last three columns are in the order of child, mother, and father.

    Issue · Github
    by Sheila

    Issue Number
    722
    State
    open
    Last Updated
    Assignee
    Array
    Milestone
    Array
  • zwangzwang torontoMember

    Do I need to upload more data?

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @zwang
    Hi,

    Sorry for the late response. I need to think about this/check with the team. I will get back to you soon.

    -Sheila

  • zwangzwang torontoMember

    Any update on this issue?

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin
    edited March 2016

    @zwang
    Hi,

    Sorry, this dropped off my radar. Thanks for the reminder.

    It seems the mom's GQ was not super high originally, and the CalculateGenotypePosteriors math showed evidence for the change. There is really strong evidence in the offspring for the variant, and the mother's GQ is relatively low. If you are concerned about this behavior, you can try lowering the de novo prior.

    -Sheila

Sign In or Register to comment.