chopsticks-vignette.R
In chopsticks: The 'snp.matrix' and 'X.snp.matrix' Classes

### R code from vignette source 'chopsticks-vignette.Rnw'

###################################################
### code chunk number 1: init
###################################################
library(chopsticks)
library(hexbin)
data(for.exercise)


###################################################
### code chunk number 2: chopsticks-vignette.Rnw:98-99
###################################################
show(snps.10)


###################################################
### code chunk number 3: chopsticks-vignette.Rnw:104-105
###################################################
summary(snp.support)


###################################################
### code chunk number 4: chopsticks-vignette.Rnw:116-117
###################################################
summary(subject.support)


###################################################
### code chunk number 5: chopsticks-vignette.Rnw:137-139
###################################################
snpsum <- summary(snps.10)
summary(snpsum)


###################################################
### code chunk number 6: plot-snpsum
###################################################
par(mfrow = c(1, 2))
hist(snpsum$MAF)
hist(snpsum$z.HWE)


###################################################
### code chunk number 7: sample-qc
###################################################
sample.qc <- row.summary(snps.10)
summary(sample.qc)


###################################################
### code chunk number 8: plot-outliners-qc
###################################################
par(mfrow = c(1, 1))
plot(sample.qc)


###################################################
### code chunk number 9: outliers
###################################################
use <- sample.qc$Heterozygosity>0
snps.10 <- snps.10[use, ]
subject.support <- subject.support[use, ]


###################################################
### code chunk number 10: if-case-control
###################################################
if.case <- subject.support$cc == 1
if.control <- subject.support$cc == 0


###################################################
### code chunk number 11: sum-case-control
###################################################
sum.cases <- summary(snps.10[if.case, ])
sum.controls <- summary(snps.10[if.control, ])


###################################################
### code chunk number 12: plot-summaries
###################################################
hb <- hexbin(sum.controls$Call.rate, sum.cases$Call.rate, xbin=50)
sp <- plot(hb)
hexVP.abline(sp$plot.vp, a=0, b=1, col="black")


###################################################
### code chunk number 13: plot-freqs
###################################################
sp <- plot(hexbin(sum.controls$MAF, sum.cases$MAF, xbin=50))
hexVP.abline(sp$plot.vp, a=0, b=1, col="white") 


###################################################
### code chunk number 14: tests
###################################################
tests <- single.snp.tests(cc, data = subject.support, snp.data = snps.10)


###################################################
### code chunk number 15: sum-tests
###################################################
summary(tests)


###################################################
### code chunk number 16: use
###################################################
use <- snpsum$MAF > 0.01 & snpsum$z.HWE^2 < 200


###################################################
### code chunk number 17: sum-use
###################################################
sum(use)


###################################################
### code chunk number 18: subset-tests
###################################################
tests <- tests[use, ]
position <- snp.support[use, "position"]


###################################################
### code chunk number 19: plot-tests
###################################################
p1 <- pchisq(tests$chi2.1df, df = 1, lower.tail = FALSE)
plot(hexbin(position, -log10(p1), xbin=50))


###################################################
### code chunk number 20: qqplot
###################################################
qq.chisq(tests$chi2.1df, df = 1)


###################################################
### code chunk number 21: more-tests
###################################################
tests <- single.snp.tests(cc, stratum, data = subject.support,
     snp.data = snps.10)
tests <- tests[use, ]
p1 <- pchisq(tests$chi2.1df, df = 1, lower.tail = FALSE)
plot(hexbin(position, -log10(p1), xbin=50))


###################################################
### code chunk number 22: more-tests-qq
###################################################
qq.chisq(tests$chi2.1df, df = 1)


###################################################
### code chunk number 23: ord
###################################################
ord <- order(p1)
top10 <- ord[1:10]
top10


###################################################
### code chunk number 24: top-10
###################################################
names <- rownames(tests)
p1[top10]
names[top10]
position[top10]


###################################################
### code chunk number 25: top10-local
###################################################
posord <- order(position)
position <- position[posord]
names <- names[posord]
local <- names[position > 9.6e+07 & position < 9.8e+07]


###################################################
### code chunk number 26: ld-2mb
###################################################
snps.2mb <- snps.10[, local]
ld.2mb <- ld.snp(snps.2mb)


###################################################
### code chunk number 27: top1
###################################################
attach(subject.support)
top <- snps.10[, "rs870041"]
top <- as.numeric(top)
glm(cc ~ top + stratum, family = "binomial")


###################################################
### code chunk number 28: top2-3
###################################################
top2 <- snps.10[, "rs7088765"]
top2 <- as.numeric(top2)
glm(cc ~ top2 + stratum, family = "binomial")
top3 <- snps.10[, "rs1916572"]
top3 <- as(top3, "numeric")
glm(cc ~ top3 + stratum, family = binomial)


###################################################
### code chunk number 29: blocks
###################################################
blocks <- split(posord, cut(position, seq(100000, 135300000, 20000)))
bsize <- sapply(blocks, length)
table(bsize)
blocks <- blocks[bsize>0]  


###################################################
### code chunk number 30: twentyfive
###################################################
posord[1:20]
blocks[1:5]


###################################################
### code chunk number 31: blocks-use
###################################################
snps.use <- snps.10[, use]
remove(snps.10)


###################################################
### code chunk number 32: mtests
###################################################
mtests <- snp.rhs.tests(cc ~ stratum, family = "binomial", 
     data = subject.support, snp.data = snps.use, tests = blocks)
summary(mtests)


###################################################
### code chunk number 33: plot-mtests
###################################################
pm <- pchisq(mtests$Chi.squared, mtests$Df, lower.tail = FALSE)
plot(hexbin(-log10(pm), xbin=50))


###################################################
### code chunk number 34: qqplot-mtests
###################################################
qq.chisq(-2 * log(pm), df = 2)