Sun, Ying [BSD] - HGD
2008-Aug-22 03:57 UTC
[R] help needed for HWE.exact in library "genetics"
Hi, I have a genotype data for both case and controls and would like to calculate the HW p-value. However, since the number of one genotype is 0, I got wired result. Would someone help me to figure it out? Or confirm it's right? Thanks a lot. ===========> library( "genetics" ) NOTE: THIS PACKAGE IS NOW OBSOLETE. The R-Genetics project has developed an set of enhanced genetics packages to replace 'genetics'. Please visit the project homepage at http://rgenetics.org for informtion.> case.data <- c("G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","A/G","A/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","A/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","A/G","G/G","A/G","G/G","G/G ","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","A/G","G/G","G/G","G/G","A/G","G/G","G/G","G/G","A/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","A/G","G/G","G/G", "G/G","G/G","G/G")> g1 <- genotype(case.data) > g1[1] "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" "G/A" "G/A" [13] "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" [25] "G/G" "G/G" "G/G" "G/G" "G/A" "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" [37] "G/G" "G/G" "G/G" "G/A" "G/G" "G/A" "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" [49] "G/G" "G/G" "G/G" "G/G" "G/A" "G/G" "G/G" "G/G" "G/A" "G/G" "G/G" "G/G" [61] "G/A" "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" [73] "G/G" "G/G" "G/A" "G/G" "G/G" "G/G" "G/G" "G/G" Alleles: G A> HWE.exact(g1)Exact Test for Hardy-Weinberg Equilibrium data: g1 N11 = 71, N12 = 9, N22 = 0, N1 = 151, N2 = 9, p-value = 1> control.data <-c("G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","A/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","A/G" ,"G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G"," G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","A/G","G/G","G/ G","G/G","G/G","G/G","G/G","A/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","A/G","G/G","G/G" ,"G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G","G/G"," G/G","G/G","G/G","G/G","G+ /G","G/G","G/G","G/G","G/G","G/G")> g2 <- genotype(control.data) > g2[1] "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" [13] "G/G" "G/G" "G/G" "G/A" "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" [25] "G/G" "G/G" "G/G" "G/G" "G/G" "G/A" "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" [37] "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" [49] "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" [61] "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" [73] "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" [85] "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" "G/A" "G/G" [97] "G/G" "G/G" "G/G" "G/G" "G/G" "G/A" "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" [109] "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" [121] "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" "G/A" "G/G" "G/G" "G/G" "G/G" [133] "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" [145] "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" [157] "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" "G/G" [169] "G/G" "G/G" "G/G" "G/G" "G/G" Alleles: G A> HWE.exact(g2)Exact Test for Hardy-Weinberg Equilibrium data: g2 N11 = 168, N12 = 5, N22 = 0, N1 = 341, N2 = 5, p-value = 1 =========== Becky This email is intended only for the use of the individua...{{dropped:12}}
You could follow the advice given when loading the library (see the code you posted for details) and use the enhanced genetics packages to cross-validate the results (ideally you should get the same answer). The results not weird though. Your working with SNPs and having a homozygote with a frequency of zero isn't that unlikely, it depends on the minor allele frequency (MAF). In this instance its the controls you appear to be concerned about, and for this locus the MAF of the minor allele is 0.0145, so in a sample of 173 individuals you would expect to see 0.0363 individuals with the recessive genotype (basic Mendelian genetics of 173 * 0.0145^2). Thus at best you might expect to see one person in a sample of that size, but its not surprising that you've not seen any. You may also find the following references of interest... Guangyong Zou, Allan Donner The Merits of Testing Hardy-Weinberg Equilibrium in the Analysis of Unmatched Case-Control Data: A Cautionary Note Annals of Human Genetics Volume 70, Issue 6, Pages 923-933 (basically advocates the use of a robust test for association such as the trend test over the two-stage design of testing for HWeqm and then testing for association). There is follow up discussion of the article in Yik Y. Teo, Andrew E. Fry, Taane G. Clark, E. S. Tai, Mark Seielstad On the Usage of HWE for Identifying Genotyping Errors Annals of Human Genetics Volume 71 Issue 5 , Pages 701 - 703 ...and a rejoinder from Zou and Donner... Guangyong Zou, Allan Donner The Reply of Zou & Donner to Teo's Commentary on their Manuscript (p 704-704) Annals of Human Genetics Volume 71 Issue 5 , Pages 701 - 703 Neil Sun, Ying [BSD] - HGD wrote:> > ===========>> library( "genetics" ) > NOTE: THIS PACKAGE IS NOW OBSOLETE. > The R-Genetics project has developed an set of enhanced genetics > packages to replace 'genetics'. Please visit the project homepage > at http://rgenetics.org for informtion. >-- View this message in context: http://www.nabble.com/help-needed-for-HWE.exact-in-library-%22genetics%22-tp19100974p19105112.html Sent from the R help mailing list archive at Nabble.com.
"Sun, Ying" asked:> I have a genotype data for both case and controls and would > like to calculate the HW p-value. However, since the > number of one genotype is 0, I got wierd result. Would > someone help me to figure it out? Or confirm it's right? > Thanks a lot....> HWE.exact(g1) > Exact Test for Hardy-Weinberg Equilibrium > data: g1 > N11 = 71, N12 = 9, N22 = 0, N1 = 151, N2 = 9, > p-value = 1Yes, that is correct. Double check it by calculating the goodness-of-fit chi-square for the same hypothesis: p <- 151/160; q <- 1-p; 80*c(p^2,2*p*q,q^2) [1] 71.253125 8.493750 0.253125 ... David Duffy. -- | David Duffy (MBBS PhD) ,-_|\ | email: davidD at qimr.edu.au ph: INT+61+7+3362-0217 fax: -0101 / * | Epidemiology Unit, Queensland Institute of Medical Research \_,-._/ | 300 Herston Rd, Brisbane, Queensland 4029, Australia GPG 4D0B994A v