R语言-计算Hardy-Weinberg Equilibrium

参考书Bioinformatics with R Cookbook

安装R包

rm(list = ls())
install.packages("GWASExactHW")
library(GWASExactHW)

生成snp数据

pA <- runif(1)
pAA <- pA^2
pAa <- 2*pA*(1-pA)
paa <- (1-pA)^2
myCounts <- rmultinom(100, 500, c(pAA, pAa, paa))
dim(myCounts)
summary(myCounts)

计算

genotypes <- data.frame(t(myCounts))
colnames(genotypes)= c("nAA", "nAa", "naa")
hwPvalues <- HWExact(genotypes)

画图

myTest <- HWExact(genotypes)
names(myTest) <- rownames(genotypes)
plot(-log10(myTest), type="b", ylab="-log10(p-value)", main="HWE p-values for SNPS")
abline(h=-log10(0.05), col="red")

输出结果

sum(myTest<0.05)
names(myTest)[which(myTest<0.05)]

P(HWE)>0.05表示群体遗传达到平衡，表示样本具有代表性，分析结论有意义