HardyWeinbergExampleSession.R
In HardyWeinberg: Statistical Tests and Graphics for Hardy-Weinberg Equilibrium

## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
knitr::opts_chunk$set(fig.align = "center")
knitr::opts_chunk$set(warnings = FALSE)
knitr::opts_chunk$set(message = FALSE)
knitr::opts_chunk$set(fig.width = 6, fig.height = 6) 

## ----preinstall---------------------------------------------------------------
#install.packages("HardyWeinberg")
library(HardyWeinberg)

## -----------------------------------------------------------------------------
x <- c(MM = 298, MN = 489, NN = 213)

## -----------------------------------------------------------------------------
maf(x)

## -----------------------------------------------------------------------------
maf(x,2)
maf(x,3)

## -----------------------------------------------------------------------------
HWTernaryPlot(x,region=0,hwcurve=FALSE,grid=TRUE,markercol="blue")

## -----------------------------------------------------------------------------
HW.test <- HWChisq(x, verbose = TRUE)

## -----------------------------------------------------------------------------
HW.test <- HWChisq(x, cc = 0, verbose = TRUE)

## -----------------------------------------------------------------------------
HW.lrtest <- HWLratio(x, verbose = TRUE)

## -----------------------------------------------------------------------------
HW.exacttest <- HWExact(x, verbose = TRUE, pvaluetype = "midp")

## -----------------------------------------------------------------------------
set.seed(123)
#HW.permutationtest <- HWPerm(x, verbose = TRUE)
#Permutation test for Hardy-Weinberg equilibrium
#Observed statistic: 0.2214896   17000 permutations. p-value: 0.6508235 

## ----alltests-----------------------------------------------------------------
#HW.results <- HWAlltests(x, verbose = TRUE, include.permutation.test = TRUE)
#                                            Statistic   p-value
#Chi-square test:                            0.2214896 0.6379073
#Chi-square test with continuity correction: 0.1789563 0.6722717
#Likelihood-ratio test:                      0.2214663 0.6379250
#Exact test with selome p-value:                    NA 0.6556635
#Exact test with dost p-value:                      NA 0.6723356
#Exact test with mid p-value:                       NA 0.6330965
#Permutation test:                           0.2214896 0.6508235

## -----------------------------------------------------------------------------
post.dens <- HWLindley(seq(-1,1,by=0.01),x)
plot(seq(-1,1,by=0.01),post.dens,type="l",xlab=expression(alpha),
     ylab=expression(pi(alpha)))
segments(0,0,0,HWLindley(0,x),lty="dotted",col="red")

## -----------------------------------------------------------------------------
HWLindley.cri(x=x)

## -----------------------------------------------------------------------------
SNP1 <- c(A=399,B=205,AA=230,AB=314,BB=107) 

## -----------------------------------------------------------------------------
HWChisq(SNP1,cc=0,x.linked=TRUE,verbose=TRUE)

## -----------------------------------------------------------------------------
HWChisq(SNP1[3:5],cc=0)

## -----------------------------------------------------------------------------
HWLratio(SNP1,x.linked=TRUE)

## -----------------------------------------------------------------------------
HWExact(SNP1,x.linked=TRUE)

## -----------------------------------------------------------------------------
HWExact(SNP1,x.linked=TRUE,pvaluetype="midp")

## -----------------------------------------------------------------------------
HWExact(SNP1[3:5],pvaluetype="midp")

## ----permutationxlinked-------------------------------------------------------
#HWPerm(SNP1,x.linked=TRUE)
#Permutation test for Hardy-Weinberg equilibrium of an X-linked marker
#Observed statistic: 7.624175   17000 permutations. p-value: 0.02211765 

## ----alltestsxlinked----------------------------------------------------------
#HWAlltests(SNP1,x.linked=TRUE,include.permutation.test=TRUE)
#                                            Statistic    p-value
#Chi-square test:                             7.624175 0.02210200
#Chi-square test with continuity correction:  7.242011 0.02675576
#Likelihood-ratio test:                       7.693436 0.02134969
#Exact test with selome p-value:                    NA 0.02085798
#Exact test with dost p-value:                      NA         NA
#Exact test with mid p-value:                       NA 0.02082957
#Permutation test:                            7.624175 0.02211765

## -----------------------------------------------------------------------------
AFtest(SNP1)

## ----bayesianx----------------------------------------------------------------
HWPosterior(males=SNP1[1:2],females=SNP1[3:5],x.linked=TRUE)

## -----------------------------------------------------------------------------
data(Markers)
Markers[1:12,]

## -----------------------------------------------------------------------------
Xt <- table(Markers[,1])
Xv <- as.vector(Xt)
names(Xv) <- names(Xt)
Xv
HW.test <- HWChisq(Xv,cc=0,verbose=TRUE)

## -----------------------------------------------------------------------------
set.seed(123)
Results <- HWMissing(Markers[,1], m = 50, method = "sample", verbose=TRUE)

## -----------------------------------------------------------------------------
set.seed(123)
Results <- HWMissing(Markers[, 1:5], m = 50, verbose = TRUE)

## -----------------------------------------------------------------------------
set.seed(123)
Results <- HWMissing(Markers[, 1:5], m = 50, statistic = "exact", verbose = TRUE)

## -----------------------------------------------------------------------------
x <- c(MM = 298, MN = 489, NN = 213)
n <- sum(x)
nM <- mac(x) 
pw4 <- HWPower(n, nM, alpha = 0.05, test = "exact", theta = 4, 
               pvaluetype = "selome")
print(pw4)
pw8 <- HWPower(n, nM, alpha = 0.05, test = "exact", theta = 8, 
               pvaluetype = "selome")
print(pw8)

## ----glyoxalasa---------------------------------------------------------------
data("Glyoxalase")
Glyoxalase <- as.matrix(Glyoxalase)
HWStrata(Glyoxalase)

## -----------------------------------------------------------------------------
pvalues <- HWExactStats(Glyoxalase)
pvalues

## -----------------------------------------------------------------------------
g1 <- c(AA=0.034, AB=0.330, BB=0.636)
g2 <- c(AA=0.349, AB=0.493, BB=0.158)
x  <- c(AA=0.270, AB=0.453, BB=0.277)

## -----------------------------------------------------------------------------
G <- cbind(g1,g2)
contributions <- HWEM(x,G=G)
contributions

## -----------------------------------------------------------------------------
p <- c(af(g1),af(g2))
contributions <- HWEM(x,p=p)
contributions

## -----------------------------------------------------------------------------
data(JPTsnps)
JPTsnps[1,]
Results <- HWPosterior(males=JPTsnps[1,1:3],
                       females=JPTsnps[1,4:6],
                       x.linked=FALSE,precision=0.05)

## -----------------------------------------------------------------------------
data(JPTsnps)
AICs <- HWAIC(JPTsnps[1,1:3],JPTsnps[1,4:6])
AICs

## ----datasets-----------------------------------------------------------------
data("CEUchr22")
CEUchr22[1:5,1:5]

## -----------------------------------------------------------------------------
Z <- MakeCounts(CEUchr22)
Z <- Z[,1:3]
head(Z)

## -----------------------------------------------------------------------------
HWTernaryPlot(Z[,1:3],patternramp=TRUE,region=0)

## -----------------------------------------------------------------------------
pminor <- maf(Z)
HWTernaryPlot(Z[pminor>0.05,],patternramp=TRUE,region=0)

## -----------------------------------------------------------------------------
hist(pminor[pminor > 0],freq=FALSE,xlab="MAF",main="CEU MAF CHR 22")

## -----------------------------------------------------------------------------
cminor <- maf(Z[pminor > 0,],option=3)
barplot(table(cminor[,1]),cex.names = 0.75)

## ----manytests----------------------------------------------------------------
Zpoly <- Z[!is.mono(Z),]
npoly <- nrow(Zpoly)
chisq.pvalues <- HWChisqStats(Zpoly,pvalues=TRUE)
exact.pvalues <- HWExactStats(Zpoly,midp=TRUE)
bonferronithreshold <- 0.05/npoly
sum(chisq.pvalues < bonferronithreshold)
sum(exact.pvalues < bonferronithreshold)

## ----ternaryCEU---------------------------------------------------------------
Zu <- UniqueGenotypeCounts(Z)
Zu <- Zu[,1:3]
HWTernaryPlot(Zu,region=1,pch=1)

## -----------------------------------------------------------------------------
HWTernaryPlot(Zu[,1:3],region=7,pch=1)

## -----------------------------------------------------------------------------
data("CEUchr22")

Z <- MakeCounts(CEUchr22)
Z <- Z[,1:3]
Z <- Z[!is.mono(Z),]

alfreq  <- af(Z)
alcount <- maf(Z,option=3)[,1]

chisq.pvals <- HWChisqStats(Z,pvalues=TRUE)

set.seed(123)
Z.sim.chi <- HWData(nm=nrow(Z),n=99,p=alfreq)
Z.sim.chi <- Z.sim.chi[!is.mono(Z.sim.chi),]
chisq.pvals.sim <- HWChisqStats(Z.sim.chi,pvalues=TRUE)

set.seed(123)
Z.sim.exa <- HWData(nm=nrow(Z),n=99,nA=alcount,conditional = TRUE)
Z.sim.exa <- Z.sim.exa[!is.mono(Z.sim.exa),]

opar <- par(mfrow=c(2,2),mar=c(3,3,2,0)+0.5,mgp=c(2,1,0))
par(mfg=c(1,1))
 qqunif(chisq.pvals,logplot = TRUE,
       plotline = 1,main="A: observed QQ uniform")
par(mfg=c(1,2))
 qqunif(chisq.pvals.sim,logplot = TRUE,xylim=15,
        main="B: simulated QQ uniform")
par(mfg=c(2,1))
 set.seed(123)
 HWQqplot(Z,nsim=5,logplot=TRUE,main="C: observed QQ HWE null")
par(mfg=c(2,2))
 set.seed(123)
 HWQqplot(Z.sim.exa,nsim=5,logplot=TRUE,main="D: simulated QQ HWE null")
par(opar)

## ----simulate-----------------------------------------------------------------
set.seed(123)
n <- 100
m <- 100
X1 <- HWData(m, n, p = rep(0.5, m))
X2 <- HWData(m, n)
X3 <- HWData(m, n, p = rep(0.5, m), f = rep(0.5, m))
X4 <- HWData(m, n, f = rep(0.5, m))
X5 <- HWData(m, n, p = rep(c(0.2, 0.4, 0.6, 0.8), 25), conditional = TRUE)
X6 <- HWData(m, n, exactequilibrium = TRUE)
opar <- par(mfrow = c(3, 2),mar = c(1, 0, 3, 0) + 0.1)
par(mfg = c(1, 1))
HWTernaryPlot(X1, main = "(a)")
par(mfg = c(1, 2))
HWTernaryPlot(X2, main = "(b)")
par(mfg = c(2, 1))
HWTernaryPlot(X3, main = "(c)")
par(mfg = c(2, 2))
HWTernaryPlot(X4, main = "(d)")
par(mfg = c(3, 1))
HWTernaryPlot(X5, main = "(e)")
par(mfg = c(3, 2))
HWTernaryPlot(X6, main = "(f)")
par(opar)

## -----------------------------------------------------------------------------
X <- HWData(nm=100,shape1=1,shape2=10)

## -----------------------------------------------------------------------------
x <- c(fA=182,fB=60,nAB=17,nOO=176)
al.fre <- HWABO(x)

## -----------------------------------------------------------------------------
x <- c(AA=20,AB=31,AC=26,BB=15,BC=12,CC=0)
#results <- HWTriExact(x)
#Tri-allelic Exact test for HWE (autosomal).
#Allele counts: A = 38 B = 73 C = 97 
#sum probabilities all outcomes 1 
#probability of the sample 0.0001122091 
#p-value =  0.03370688 

## -----------------------------------------------------------------------------
x <- c(AA=20,AB=31,AC=26,BB=15,BC=12,CC=0)
x <- toTriangular(x)
x
m <- c(A=0,B=0,C=0)
results <- HWNetwork(ma=m,fe=x)

## -----------------------------------------------------------------------------
set.seed(123)
x <- c(AA=20,AB=31,AC=26,BB=15,BC=12,CC=0)
x <- toTriangular(x)
#results <- HWPerm.mult(x)
#Permutation test for Hardy-Weinberg equilibrium (autosomal).
#3 alleles detected.
#Observed statistic: 0.0001122091   17000 permutations. p-value: 0.03405882 

## -----------------------------------------------------------------------------
males   <- c(A=1,B=21,C=34) 
females <- c(AA=0,AB=1,AC=0,BB=8,BC=24,CC=15)
results <- HWTriExact(females,males)

## -----------------------------------------------------------------------------
males   <- c(A=1,B=21,C=34) 
females <- toTriangular(c(AA=0,AB=1,AC=0,BB=8,BC=24,CC=15))
results <- HWNetwork(ma=males,fe=females)

## -----------------------------------------------------------------------------
data("NistSTRs")
NistSTRs[1:5,1:5]
n <- nrow(NistSTRs)
p <- ncol(NistSTRs)/2
n
p

A1 <- NistSTRs[,1]
A2 <- NistSTRs[,2]
GenotypeCounts <- AllelesToTriangular(A1,A2)
print(GenotypeCounts)

set.seed(123)
#out <- HWPerm.mult(GenotypeCounts)
#Permutation test for Hardy-Weinberg equilibrium (autosomal).
#7 alleles detected.
#Observed statistic: 2.290724e-11   17000 permutations. p-value: 0.8644706 

## -----------------------------------------------------------------------------
#Results <- HWStr(NistSTRs,test="permutation")
#29 STRs detected.
#> Results
#          STR   N Nt MinorAF MajorAF     Ho     He     Hp   pval
#1    CSF1PO-1 361  7  0.0055  0.3601 0.7341 0.7194 1.4016 0.8614
#2  D10S1248-1 361  9  0.0014  0.3075 0.7645 0.7586 1.5552 0.6488
#3   D12S391-1 361 16  0.0014  0.1717 0.8975 0.8909 2.3686 0.8955
#4   D13S317-1 361  8  0.0014  0.3255 0.7618 0.7837 1.7102 0.1049
#5   D16S539-1 361  7  0.0180  0.3144 0.7645 0.7600 1.5933 0.1641
#6    D18S51-1 361 15  0.0014  0.1704 0.8560 0.8758 2.2139 0.1926
#7   D19S433-1 361 15  0.0014  0.3615 0.7673 0.7694 1.7624 0.9667
#8   D1S1656-1 361 15  0.0028  0.1496 0.9252 0.8992 2.4073 0.8699
#9    D21S11-1 361 16  0.0014  0.2825 0.8227 0.8288 1.9946 0.6156
#10 D22S1045-1 361  8  0.0055  0.3823 0.7535 0.7220 1.4823 0.9785
#11  D2S1338-1 361 12  0.0014  0.1856 0.8698 0.8814 2.2464 0.5066
#12   D2S441-1 361 11  0.0014  0.3435 0.7867 0.7693 1.6734 0.7125
#13  D3S1358-1 361  9  0.0014  0.2729 0.7562 0.7900 1.6437 0.3789
#14   D5S818-1 361  9  0.0014  0.3878 0.7064 0.6977 1.3939 0.7615
#15  D6S1043-1 361 14  0.0014  0.2964 0.7978 0.8232 2.0059 0.0220
#16   D7S820-1 361  9  0.0014  0.2562 0.8310 0.8161 1.7925 0.5045
#17  D8S1179-1 361 10  0.0014  0.3296 0.7839 0.8072 1.8386 0.8606
#18   F13A01-1 361 12  0.0014  0.3490 0.7590 0.7353 1.5471 0.8793
#19     F13B-1 361  6  0.0055  0.3892 0.7008 0.7175 1.3790 0.8276
#20   FESFPS-1 361  7  0.0014  0.4114 0.6731 0.6930 1.3085 0.8506
#21      FGA-1 361 14  0.0014  0.2050 0.8670 0.8594 2.1163 0.0528
#22      LPL-1 361  8  0.0014  0.4224 0.7175 0.6947 1.3386 0.9721
#23  Penta_C-1 361 10  0.0014  0.3947 0.7701 0.7523 1.6035 0.0248
#24  Penta_D-1 361 13  0.0014  0.2327 0.8587 0.8247 1.8925 0.1464
#25  Penta_E-1 361 19  0.0014  0.1994 0.8920 0.8910 2.4391 0.4811
#26     SE33-1 361 39  0.0014  0.0942 0.9501 0.9483 3.1472 0.2506
#27     TH01-1 361  8  0.0014  0.3449 0.7424 0.7646 1.5616 0.2379
#28     TPOX-1 361  8  0.0014  0.5249 0.6537 0.6404 1.2572 0.5545
#29      vWA-1 361 10  0.0014  0.2839 0.8061 0.8076 1.7577 0.4501