Forum Login Issue:
Currently the "Log in with Google" button redirects you to a "Page not found." This is an issue that our forum vendors are working on fixing. In the meantime, while on the "Page not found" you can edit the URL to delete the second gatk, firecloud, or wdl (depending on what subforum you are acessing).
ex: https://gatkforums.broadinstitute.org/gatk/gatk/entry/...

--input-prior default value?

Issue · Github
by Sheila

Issue Number
1595
State
closed
Last Updated
Assignee
Array
Milestone
Array
Closed By
vdauwera

Best Answers

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie
    Accepted Answer

    Alright, never mind anything I said above -- after a bit of digging it's a bit more complicated than I thought, and not a version/model related thing. The priors you quote from the article are derived from the SNP (or indel) heterozygosity (each a scalar) using a pre-determined equation. If you don't want to follow the pre-determined equation, then you use the input_prior argument, which according to the dev I got this from should override the heterozygosity argument. There's some notes about how heterozygosity is used here.

Answers

  • tommycarstensentommycarstensen United KingdomMember
    edited January 2017

    Thanks a lot for this @Geraldine_VdAuwera! It makes sense. Your current documentation suggests to the reader, that the priors are not flat:

    For completely flat priors, specify the same value (=1/(2N+1)) 2N times, e.g. -inputPrior 0.33 -inputPrior 0.33 for the single-sample diploid case.

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie

    Ho hmm that's a fair point. I have a Q30 confidence that they are in fact flat, but will check and clarify the doc.

  • tommycarstensentommycarstensen United KingdomMember
    edited January 2017

    Q30 :) Not quite sure what that translates to on the VdA confidence scale :)

    I'm violating the copyright laws by pasting here from the supplementary to the 2016 Nature paper "The Simons Genome Diversity Project: 300 genomes from 142 diverse populations" by Mallick, et.al., which seems to suggest the priors are not flat by default:

    http://www.nature.com/nature/journal/v538/n7624/extref/nature18964-s1.pdf

    Most analyses in this paper are based on single-sample genotypes determined using a reference-bias free modification of GATK. We did not perform multi-sample genotyping as we were concerned that this could induce biases in population genetic analyses. Specifically, we were concerned that the GATK UnifiedGenotyper has a built-in prior for Bayesian SNP calling that assumes that the site is more likely to be homozygous for the reference allele than homozygous for the variant allele. For a diploid sample, the default priors for a homozygous reference, heterozygote and homozygous non-reference genotypes are (0.9985, 0.001, 0.0005), respectively. When there is ambiguity in a heterozygote, GATK prefers the reference homozygote. This is a reference bias, and while this bias is not typically problematic for medical studies, it can complicate interpretation of population genetics signals. With the Genome Sequencing and Analysis Group at the Broad Institute, we developed an alternative model that was integrated into the UnifiedGenotyper, allowing reference-bias free priors to be specified. We are using a prior (0.4995, 0.001, 0.4995). Details are at: https://www.broadinstitute.org/gatk/gatkdocs/org_broadinstitute_gatk_tools_walkers_genotyper_UnifiedGenotyper.php#--input_prior.

  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie
    Hah. Don't put too much stock in it, it's not a well-calibrated score ;)

    Well this sounds like it might be a bit dated. We haven't been called GSA for over two years;
    that was back before we were hybridized with other teams including Picard devs to change the world (and I promised you butterflies and the incubation has taken massively longer than I thought but we're almost there and that's a story for another time). The genotyping model has gone through several iterations since I've been with the team (the code is littered with their carcasses -- I mean models, not team members); and certainly the population genetics thing has become very important to us (hello exac). Do they say what version they used? Kind of unfortunate to link to a web doc that gets updated with every version.

    Ah you helpfully linked to the supplementary materials, allowing me to see they used version 2.5. Well that's not recent, blimey. Can't say for certain yet (will have to go to the code for that) but that's shaping up to be a likely explanation. The passage of time and the rotation of astral bodies...
  • Geraldine_VdAuweraGeraldine_VdAuwera Cambridge, MAMember, Administrator, Broadie
    Accepted Answer

    Alright, never mind anything I said above -- after a bit of digging it's a bit more complicated than I thought, and not a version/model related thing. The priors you quote from the article are derived from the SNP (or indel) heterozygosity (each a scalar) using a pre-determined equation. If you don't want to follow the pre-determined equation, then you use the input_prior argument, which according to the dev I got this from should override the heterozygosity argument. There's some notes about how heterozygosity is used here.

Sign In or Register to comment.