Attention:
The frontline support team will be slow on the forum because we are occupied with the GATK Workshop on March 21st and 22nd 2019. We will be back and more available to answer questions on the forum on March 25th 2019.

Unexpected genotype likelihoods & genotype qualities

I noted that some of my genotype qualities are equal to 0. I interpret this as if two genotypes were equally likely. However, sometimes this interpretation does not look consistent with the primary data. For instance, I am confused with the example below:

GT:AD:DP:GQ:PL 0/0:5,0:5:15:0,15,171
GT:AD:DP:GQ:PL 0/0:21,0:21:0:0,0,430

Both genotypes are homozygous REFs.

The first one is supported by 5 out of 5 reads. It scores genotype quality 15. This looks reasonable: there is some probability of heterozygous genotype because of the low coverage.

The second genotype is supported by 21 of 21 reads. However, it suddenly yields the genotype score of … zero ? How could it be that the heterozygous is as likely as homozygous given such evidence?

Investigating this example further, I noted that there was a significant difference in ALT allele frequency between these two variants in the whole studied dataset (~500 samples). The first variant has quite low ALT frequency (<1%). The second variant has quite frequent ALT allele ( ~50%). Can this explain the unexpected zero genotype quality calculated for the second genotype?

What factors are taken into account when calculating the genotype likelihoods?
Is this example an expected output?
Could it be an error?

Thank you in advance,

Best Answers

Answers

  • Dear Sheila

    Thank you for the prompt reply.
    Please see below the requested information about my commands and gatk version.

    I have an example of vcf file, which contains quite a number of zero GQ cases.
    The ftp link is at the end of this message, I may keep it alive for some days.
    Please let me know when/if you do not need it, so I will disable the link.
    Sorry, the file is quite large.

    In practical terms I would not like having genotypes ./. in my examples above...
    I would rather prefer to understand how GQ could be 0 when 21 reads of 21 support 0/0?
    Or where I was wrong in generating or interpreting the data?

    Any other comments on my commands and data would also be appretiated.

    Many thanks again,
    Best regards,
    Alexey

    gatk version

    gatk-3.4-46

    making gvcfs

    "${java7}" -Xmx60g -jar "${gatk}" \
    -T HaplotypeCaller \
    -R "${ref_genome}" \
    -L "${nextera_targets_intervals}" -ip 100 \
    -I "${idr_bqr_bam}" \
    -o "${gvcf}" \
    -ERC GVCF \
    -nct 12 \
    2> "${gvcf_log}"

    combining gvcfs

    "${java7}" -Xmx60g -jar "${gatk}" \
    -T CombineGVCFs \
    -R "${ref_genome}" \
    -L "${nextera_targets_intervals}" -ip 100 \
    -V "${source_gvcfs}" \
    -o "${combined_gvcf}" \
    2> "${combining_gvcf_log}"

    genotyping combined gvcfs

    "${java7}" -Xmx60g -jar "${gatk}" \
    -T GenotypeGVCFs \
    -R "${ref_genome}" \
    -L "${nextera_targets_intervals}" -ip 100 \
    -maxAltAlleles 25 \
    -V "${source_gvcfs}" \
    -o "${raw_vcf}" \
    -nt 14 \
    &> "${GenotypeGVCFs_log}"

    Example of data: vcf file of ~7.8Gb

    ftp host: mgqnap.medschl.cam.ac.uk
    user: sheila
    password: sheila
    I will keep access alive till Sunday

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    @alexey_larionov, we have a specific method for accepting test data for bug reports. We cannot take full files, so you will need to subset the data to a small region around the site of interest. We will need you to provide original read data from the bam file as the VCFs are not sufficient for us to investigate this problem. Also, we need you to upload your data to our FTP. See this FAQ article for full details.

  • Dear Geraldine,
    Thank you for your reply,

    I will read the instructions and try to follow the procedure within a week.

    Best regards,
    Alexey

    PS: I have disabled the ftp link

  • laehnemannlaehnemann University of DuesseldorfMember

    Hi there,

    I can only second Alexey's question. I have the same trouble understanding the double zero PL values when there is (very) clear evidence for only one genotype. Can this evidence within the sample be outweighed so heavily by evidence in other samples? Can this be avoided by only calling variants on samples with a very close relationship and thus little genotype diversity?!
    Any hints and pointers would be welcome!

    Best regards,
    David

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    @laehnemann No, this is not an expected effect of joint analysis. It may be a bug (hence our request for test data) or an effect of realignment. Have you tried generating the output bam to see what HaplotypeCaller sees? Also, which version are you using?

  • laehnemannlaehnemann University of DuesseldorfMember

    I can go and generate the output bam, if that helps. I can also go and follow the bug report guide, if you'd like.

    I used version 3.4-46-gbc02625 of the GATK for both HaplotypeCaller and ApplyRecalibration.

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    Yes please, I don't think we ever received test data from Alexey.

  • laehnemannlaehnemann University of DuesseldorfMember

    OK, I'll identify some spots and compile some testing data tomorrow!

  • Thank you for all who took time to look at this issue for me. I think that I now have understood the initial Sheila’s explanations and I fully accept the fix (emitting “./.” in zero-quality genotype cases).

    This is my current version of what happened before the fix:

    My initial concern was not about equal likelihoods in two genotypes. My main concern was that the likelihoods were not consistent with the alleles depth. It looked like a possible bug affecting ~1% of genotypes in my data. The thread provided by asilver explained that this was not a bug, but a way of reporting the checks, which HC makes after calculating likelihoods. It seems that the two most common scenarios were

    1) In case if there is no evidence for ALT allele at the site (i.e. the raw likelihood for HOM-REF is the highest) HC checks whether there are other, non-ALT, mismatches at the site. If HC finds the other mismatches it caps HOM-REF likelihood, levelling it with HET, making the HOM-REF genotype being emitted with zero quality. Importantly, this was not intended to imply any possibility of the HET. This was just a way of saying that HC could not be sure in the HOM-REF either.

    2) In case if the raw likelihood for HET is the highest (i.e. there is evidence for ALT allele at the site) HC looks whether there are some soft-clipped basses near the site. Again, if the soft-clipped bases are found, the HET likelihoods are capped to HOM-REF. Before the fix this could result to HOM-REF genotype with zero quality. Again, this did not imply presence of HOM-REF, this was a way of saying that we cannot be sure in HET either.

    The mismatches in vicinity of a variant explained all checked zero-quality genotypes in my data (I was able to check only selected few, of course).

    Again thank you everyone for taking time to explain this to me (and for stopping emitting genotypes with zero qualities :)

    PS: and I think that there is no need for me to upload any data at this point :))

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @alexey_larionov
    Hi,

    I am happy you are happy with the explanation. However, the ./. genotype has only been implemented in GVCF mode when the GQ is 0. If you have test files that recreate the issue in normal HaplotypeCaller mode (without using -ERC GVCF), we can put in the fix for that case as well.

    Thanks,
    Sheila

  • laehnemannlaehnemann University of DuesseldorfMember

    Sorry for this really late reply, but I have finally compiled a (well, not all that) little snippet of human chromosome 8 (17e6 to 20e6) to demonstrate my problem, i.e. that after joined GenotypeGVCFs on individual sample *.g.vcf files from HaplotypeCaller, I have calls where:

    1. The GQ is 0, i.e. two genotypes have the same likelihood or one of the above explained heuristics for setting GQ to zero kicks in.
    2. A GT is called, as in the GT field is not "./.".
    3. The AD is very high in at least one allele, pointing to a very clear GT with proper read support.

    The only idea I have so far comes from looking at the (HaplotypeCaller) output bam files, where according positions in according samples are not covered by any reads any more. And it is that HaplotypeCaller does the initial GT call on the read alignment from the input bam file, then bases the GQ on the new local assembly of reads and doesn't update the GT or the AD field(s) to accomodate the new read coverage -- or rather lack thereof.

    Long story short, have a look at this archive on your ftp server, all necessary inputs, outputs and logs are included (I hope):
    unexpected_0GQ_with_high_AD.tar.gz

    P.S.: I have kept the input slice bam and could easily rerun the same analysis with the newest GATK version 3.5, if there is any interest...

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie admin

    @laehnemann can you please check whether this still happens with the latest nightly build?

  • laehnemannlaehnemann University of DuesseldorfMember

    Hi again, @Geraldine_VdAuwera

    sorry for the slow response -- apparently I didn't have eMail-notifications for the forum set up to the sensitivity I should have...

    I just ran my snippet with the last nightly build and am sorry to say, but the same still happens. When filtering for sites where at least one sample has all of the following three:

    1. A genotype call.
    2. A GQ of 0.
    3. One of the alleles with a coverage of at least 30 (FMT/AD).

    I get the exact same 7 sites in the snippet of chr8 that I got before -- with the exact same information in the FORMAT fields. Only the annotations in the INFO FIELDS seem to differ:
    1. A new INFO annotation appeared: ExcessHet
    2. All RankSum INFO annotations seem to have values that differ quite drastically between the versions I tested (3.4-46-gbc02625 and 3.5-2016-03-18-gcfe332a), e.g. for the last site in the file:

    old: BaseQRankSum=2.87;ClippingRankSum=0.067;MQRankSum=-0.762;QD=34.11;ReadPosRankSum=0.762
    new: BaseQRankSum=9;ClippingRankSum=9;MQRankSum=9;QD=34.11;ReadPosRankSum=104.23
    

    And looking at the example sites in the output bam files of the respective samples with samtools tview I have seen that in one case, the bam does not have the coverage that the FMT/AD values suggest. In another case it does, and for that case mapping quality also seems fine -- so no explanation for the GQ=0 there.

    In any case, I bundled all of the new output, including a README of what I did, into unexpected_0GQ_with_high_AD-GATK_nightly-2016-03-18-gcfe332a.tar.gz, which should now be on the public ftp server. Let me know what you make of it...

    Best regards,
    David

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @laehnemann
    Hi David,

    Thanks for submitting test files. I am just testing them out now.

    -Sheila

    Issue · Github
    by Sheila

    Issue Number
    723
    State
    closed
    Last Updated
    Assignee
    Array
    Milestone
    Array
    Closed By
    chandrans
  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @laehnemann
    Hi David,

    I have not been able to find sites with AD values greater than 30 that have GQ0. Can you tell me an exact position I should look at?

    I did find some cases where the GQ is 0 and there is a genotype that is not a no-call. This only happens when I call variants using HaplotypeCaller in normal mode on all samples. However, this does not happen when I use GVCF mode (then, the samples with GQ0 are set to no-call). I am about to submit a bug report for that case.

    -Sheila

  • laehnemannlaehnemann University of DuesseldorfMember

    Hi @Sheila ,

    for the files in the uploaded in the *.tar.gz-files on your ftp server, I used bcftools 1.3 to find such sites, using the following filter command line (also see the README_STEPS file in the uploaded bundle):

    bcftools filter -i 'FMT/GT!="."&FMT/GQ=0&(FMT/AD[0]>30|FMT/AD[1]>30)' all.chr8_17e6-20e6.vcf >all.chr8_17e6-20e6.GTcalled_GQ-0_AD-30.vcf
    

    Thus, the bundle I uploaded to the ftp-server already contains a selection of such sites in the following file:
    all.chr8_17e6-20e6.GTcalled_GQ-0_AD-30.vcf

    Else, you could use the same filter command on any file produced by HaplotypeCaller on a diploid genome, to find such sites -- or increase the sensitivity by decreasing the required coverage from 30, e.g. down to 15.

    Otherwise I just realised, that the log-files in the newer bundle I uploaded (using the nightly build of GATK 3.5) are empty. Sorry for that, but in essence they should be the same as the ones from the earlier bundle (unexpected_0GQ_with_high_AD.tar.gz). In any case, I ran HaplotypeCaller in GVCF...

    Hope this helps?

    Best regards,
    David

  • SheilaSheila Broad InstituteMember, Broadie, Moderator admin

    @laehnemann
    Hi David,

    Indeed your VCF helped a lot. I will discuss with the team and get back to you.

    -Sheila

  • laehnemannlaehnemann University of DuesseldorfMember

    Hi @Sheila

    thanks for the update, I'm curious if you find something there.

    Best regards,
    David

  • laehnemannlaehnemann University of DuesseldorfMember

    @Geraldine_VdAuwera Thanks for the Update! Will make sure to upgrade to the new release, once it's there.

  • FabienneFabienne ParisMember

    Hello,

    I'm using GATK version 3.7 and I have a simlar problem.
    I've used HaplotypeCaller to generate inidvidual GVCFs and after GenotypeGVCFs on ultiple individual samples.
    In the resulting vcf, I have some genotypes 0/0 with a GQ=0, whereas I was waiting for genotypes ./. when 2 genotypes have the same probability.

    Could you please tell me how to do to force GATK to put ./. when GQ=0 ?

    Thank very much.
    Fabienne

  • bhanuGandhambhanuGandham Cambridge MAMember, Administrator, Broadie, Moderator admin

    Hi Fabienne,

    You could "SelectVariants" with --set-filtered-gt-to-nocall for this purpose.
    Also, here is a document explaining the feature: https://software.broadinstitute.org/gatk/documentation/article?id=12350

    Please let me know if this helps. :)

    Regards
    Bhanu

  • FabienneFabienne ParisMember

    Hi Bhanu,

    Thank you for your answer. Yes, it could help.
    It's a way to get around the problem of calling of GenotypeGVCFs. Do you know if this way of working is considered as a bug by the GATK team ?

    Thank you.

    Fabienne

  • bhanuGandhambhanuGandham Cambridge MAMember, Administrator, Broadie, Moderator admin

    Hi Fabienne,

    To look into this further would you please give me the exact command you are using and a few records from the resulting vcf with genotypes 0/0 with a GQ=0.
    Thank you.

    Regards
    Bhanu

Sign In or Register to comment.