Simple question pvalHigh

[Bartosz Ulaszewski]

Dear Mike and all AfteRAD users,
I have a simple question about the pvalHigh. In example shown in the manual we have two potential alleles one with 350 and second 500 which gives a p-value of 1.50-e7, which is below the threshold. Does it mean that both of these potential alleles will be treated as separate homozygous?

[Mike Sovic]

Hi Bartosz,

In this example scenario, the p-value associated with the binomial test is 1.5e-7.  If this value is lower than the relevant critical p-value (pvalHigh in this case, default 0.00001), then this sample would be scored as a homozygote for the allele with 500 reads, and the allele with 350 reads would just be ignored (it is not treated separately as a homozygote).  Clearly, this is not optimal, as this sample should be scored as a heterozygote.  

This is actually something I've been thinking a lot about just recently, trying to come up with a way to reduce this bias toward homozygous genotype calls when read depths get high.  One thought I have is to add a check for each genotype, where the locus would be scored as heterozygous if it meets either of two criteria:

1.)  the p-value criterion already used, or

2.)  a criterion in which the allele with the lower of the two counts (350 in the example above) is simply some proportion of the higher of the two.  For example, if it is at least 20% of the higher (350 is much more than 20% of 500), then the locus would be scored as heterozygous, even if it doesn't meet the p-value criterion.

I have some replicate datasets that I hope to use soon to try out a few options such as this to see if we can improve the genotyping method a bit.

Note that this really only seems to become an issue if your read depths at the locus being scored are high (i.e. in the 100's).  There are a couple ways you can check your datasets to get some sense of the level of this effect…

1.)  Run replicate samples (libraries), and check the repeatability levels of your genotypes.  There is a script available to score repeatability.

2.)  In the absence of replicate samples, one way to get some sense of this is to browse through some of the files in the folder TempFiles/RawReadCountFiles (I'd recommend specifically the ones with "NonParalogous" in the name).  There is one of these files for each sample in your dataset.  The structure of each of these files is all the same (same loci/alleles), but the counts for each allele of course apply to the specific sample.  These counts are those used for the genotyping.  So, you can compare these counts to the actual called genotypes (for example, in the Haplotypes or SNPMatrix files) to make sure you're not seeing obviously heterozygous samples (based on the read counts) incorrectly genotyped as homozygous.  If you see this happening, you might consider adjusting those critical p-values - esp. the one for pvalHigh.

If anyone has any suggestions/thoughts for ways to improve this, they are certainly welcome.

[Bartosz Ulaszewski]

Hi Mike,
Thanks for the answer, I'll try to check my data as you suggested, maybe something will come to my mind to improve this step.