Calculate HWE for multiple SNPs

hwe.R

# Create some SNPs:

px <- 0.5
py <- 0.3
pz <- 0.2

snpx <- rbinom(1000, 2, px)
snpy <- rbinom(1000, 2, py)
snpz <- rbinom(1000, 2, pz)

allelescore <- snpx + snpy + snpz


# Exact way is to perform a multinomial test, 
# and to calculate exact expected probabilities for each count
# There must be an algorithmic way to do this but I just wrote it out instead

xx <- px^2
Xx <- 2*px*(1-px)
XX <- (1-px)^2

yy <- py^2
Yy <- 2*py*(1-py)
YY <- (1-py)^2

zz <- pz^2
Zz <- 2*pz*(1-pz)
ZZ <- (1-pz)^2


g0 <- 
XX * YY * ZZ

g1 <- 
XX * YY * Zz +
XX * Yy * ZZ + 
Xx * YY * ZZ

g2 <- 
XX * YY * zz +
XX * yy * ZZ +
xx * YY * ZZ +
XX * Yy * Zz +
Xx * YY * Zz +
Xx * Yy * ZZ

g3 <- 
Xx * Yy * Zz +
XX * Yy * zz +
XX * yy * Zz +
Xx * YY * zz +
Xx * yy * ZZ +
xx * YY * Zz +
xx * Yy * ZZ

g4 <- 
XX * yy * zz +
xx * YY * zz +
xx * yy * ZZ +
Xx * Yy * zz +
Xx * yy * Zz +
xx * Yy * Zz

g5 <-
Xx * yy * zz +
xx * Yy * zz +
xx * yy * Zz

g6 <- 
xx * yy * zz


g <- c(g0, g1, g2, g3, g4, g5, g6)

# Check that sum of probabilities for each number of alleles is 1
sum(g)

g*1000
table(allelescore)
obs <- as.numeric(table(allelescore))


# Perform multinomial test
library(EMT)
multinomial.test(obs=obs, p=g, MonteCarlo=TRUE, ntrial=1e8)


## Very slow, imprecise, alternative method is to perform meta analysis of HWE test for each SNP, e.g.

library(HardyWeinberg)

tx <- HWExact(table(snpx))
ty <- HWExact(table(snpy))
tz <- HWExact(table(snpz))

Stouffer.test <- function(p, w) { # p is a vector of p-values
	if (missing(w)) {
		w <- rep(1, length(p))/length(p)
	} else {
		if (length(w) != length(p))
		stop("Length of p and w must equal!")
	}
	Zi <- qnorm(1-p)
	Z  <- sum(w*Zi)/sqrt(sum(w^2))
	p.val <- 1-pnorm(Z)
	return(c(Z = Z, p.value = p.val))
}

Stouffer.test(c(tx$pval, ty$pval, tz$pval))