###################################################
### chunk number 1: init
###################################################
#line 76 "vignettes/snpStats/inst/doc/snpStats-vignette.Rnw"
require(snpStats)
require(hexbin)
data(for.exercise)


###################################################
### chunk number 2: 
###################################################
#line 106 "vignettes/snpStats/inst/doc/snpStats-vignette.Rnw"
show(snps.10)


###################################################
### chunk number 3: 
###################################################
#line 113 "vignettes/snpStats/inst/doc/snpStats-vignette.Rnw"
summary(snp.support)


###################################################
### chunk number 4: 
###################################################
#line 125 "vignettes/snpStats/inst/doc/snpStats-vignette.Rnw"
summary(subject.support)


###################################################
### chunk number 5: 
###################################################
#line 147 "vignettes/snpStats/inst/doc/snpStats-vignette.Rnw"
summary(snps.10)


###################################################
### chunk number 6: 
###################################################
#line 153 "vignettes/snpStats/inst/doc/snpStats-vignette.Rnw"
snpsum <- col.summary(snps.10)
summary(snpsum)


###################################################
### chunk number 7: plot-snpsum
###################################################
#line 161 "vignettes/snpStats/inst/doc/snpStats-vignette.Rnw"
par(mfrow = c(1, 2))
hist(snpsum$MAF)
hist(snpsum$z.HWE)


###################################################
### chunk number 8: sample-qc
###################################################
#line 176 "vignettes/snpStats/inst/doc/snpStats-vignette.Rnw"
sample.qc <- row.summary(snps.10)
summary(sample.qc)


###################################################
### chunk number 9: plot-outliners-qc
###################################################
#line 184 "vignettes/snpStats/inst/doc/snpStats-vignette.Rnw"
par(mfrow = c(1, 1))
plot(sample.qc)


###################################################
### chunk number 10: outliers
###################################################
#line 193 "vignettes/snpStats/inst/doc/snpStats-vignette.Rnw"
use <- sample.qc$Heterozygosity>0
snps.10 <- snps.10[use, ]
subject.support <- subject.support[use, ]


###################################################
### chunk number 11: if-case-control
###################################################
#line 205 "vignettes/snpStats/inst/doc/snpStats-vignette.Rnw"
if.case <- subject.support$cc == 1
if.control <- subject.support$cc == 0


###################################################
### chunk number 12: sum-case-control
###################################################
#line 211 "vignettes/snpStats/inst/doc/snpStats-vignette.Rnw"
sum.cases <- col.summary(snps.10[if.case, ])
sum.controls <- col.summary(snps.10[if.control, ])


###################################################
### chunk number 13: plot-summaries
###################################################
#line 217 "vignettes/snpStats/inst/doc/snpStats-vignette.Rnw"
hb <- hexbin(sum.controls$Call.rate, sum.cases$Call.rate, xbin=50)
sp <- plot(hb)
hexVP.abline(sp$plot.vp, 0, 1, col="black")


###################################################
### chunk number 14: plot-freqs
###################################################
#line 226 "vignettes/snpStats/inst/doc/snpStats-vignette.Rnw"
sp <- plot(hexbin(sum.controls$MAF, sum.cases$MAF, xbin=50))
hexVP.abline(sp$plot.vp, 0, 1, col="white")


###################################################
### chunk number 15: tests
###################################################
#line 245 "vignettes/snpStats/inst/doc/snpStats-vignette.Rnw"
tests <- single.snp.tests(cc, data = subject.support, snp.data = snps.10)


###################################################
### chunk number 16: sum-tests
###################################################
#line 260 "vignettes/snpStats/inst/doc/snpStats-vignette.Rnw"
summary(tests)


###################################################
### chunk number 17: use
###################################################
#line 276 "vignettes/snpStats/inst/doc/snpStats-vignette.Rnw"
use <- snpsum$MAF > 0.01 & snpsum$z.HWE^2 < 200


###################################################
### chunk number 18: sum-use
###################################################
#line 283 "vignettes/snpStats/inst/doc/snpStats-vignette.Rnw"
sum(use)


###################################################
### chunk number 19: subset-tests
###################################################
#line 288 "vignettes/snpStats/inst/doc/snpStats-vignette.Rnw"
tests <- tests[use]
position <- snp.support[use, "position"]


###################################################
### chunk number 20: plot-tests
###################################################
#line 296 "vignettes/snpStats/inst/doc/snpStats-vignette.Rnw"
p1 <- p.value(tests, df=1)
plot(hexbin(position, -log10(p1), xbin=50))


###################################################
### chunk number 21: qqplot
###################################################
#line 303 "vignettes/snpStats/inst/doc/snpStats-vignette.Rnw"
chi2 <- chi.squared(tests, df=1)
qq.chisq(chi2,  df = 1)


###################################################
### chunk number 22: more-tests
###################################################
#line 328 "vignettes/snpStats/inst/doc/snpStats-vignette.Rnw"
tests <- single.snp.tests(cc, stratum, data = subject.support,
     snp.data = snps.10)
tests <- tests[use]
p1 <- p.value(tests, df = 1)
plot(hexbin(position, -log10(p1), xbin=50))


###################################################
### chunk number 23: more-tests-qq
###################################################
#line 335 "vignettes/snpStats/inst/doc/snpStats-vignette.Rnw"
chi2 <- chi.squared(tests, df=1)
qq.chisq(chi2, df = 1)


###################################################
### chunk number 24: ord
###################################################
#line 348 "vignettes/snpStats/inst/doc/snpStats-vignette.Rnw"
ord <- order(p1)
top10 <- ord[1:10]
top10


###################################################
### chunk number 25: top-10
###################################################
#line 357 "vignettes/snpStats/inst/doc/snpStats-vignette.Rnw"
names <- tests@snp.names
p1[top10]
names[top10]
position[top10]


###################################################
### chunk number 26: top10-local
###################################################
#line 370 "vignettes/snpStats/inst/doc/snpStats-vignette.Rnw"
posord <- order(position)
position <- position[posord]
names <- names[posord]
local <- names[position > 9.6e+07 & position < 9.8e+07]


###################################################
### chunk number 27: top1
###################################################
#line 391 "vignettes/snpStats/inst/doc/snpStats-vignette.Rnw"
cc <- subject.support$cc
stratum <- subject.support$stratum
top <- as(snps.10[, "rs870041"], "numeric")
glm(cc ~ top + stratum, family = "binomial")


###################################################
### chunk number 28: top2
###################################################
#line 400 "vignettes/snpStats/inst/doc/snpStats-vignette.Rnw"
top2 <- as(snps.10[, "rs10882596"], "numeric")
fit <- glm(cc ~ top2 + stratum, family = "binomial")
summary(fit)


###################################################
### chunk number 29: estimates
###################################################
#line 423 "vignettes/snpStats/inst/doc/snpStats-vignette.Rnw"
localest <- snp.rhs.estimates(cc~stratum, family="binomial", sets=local,
                              data=subject.support, snp.data=snps.10)


###################################################
### chunk number 30: list-estimates
###################################################
#line 431 "vignettes/snpStats/inst/doc/snpStats-vignette.Rnw"
localest[1:5]
localest["rs10882596"]


###################################################
### chunk number 31: fast-estimates
###################################################
#line 446 "vignettes/snpStats/inst/doc/snpStats-vignette.Rnw"
allest <- snp.rhs.estimates(cc~strata(stratum), family="binomial", sets=use,
                              data=subject.support, snp.data=snps.10)
length(allest)


###################################################
### chunk number 32: check-estimates
###################################################
#line 453 "vignettes/snpStats/inst/doc/snpStats-vignette.Rnw"
allest["rs10882596"]


###################################################
### chunk number 33: blocks
###################################################
#line 477 "vignettes/snpStats/inst/doc/snpStats-vignette.Rnw"
blocks <- split(posord, cut(position, seq(100000, 135300000, 20000)))
bsize <- sapply(blocks, length)
table(bsize)
blocks <- blocks[bsize>0]  


###################################################
### chunk number 34: twentyfive
###################################################
#line 486 "vignettes/snpStats/inst/doc/snpStats-vignette.Rnw"
posord[1:20]
blocks[1:5]


###################################################
### chunk number 35: blocks-use
###################################################
#line 496 "vignettes/snpStats/inst/doc/snpStats-vignette.Rnw"
snps.use <- snps.10[, use]
remove(snps.10)


###################################################
### chunk number 36: mtests
###################################################
#line 503 "vignettes/snpStats/inst/doc/snpStats-vignette.Rnw"
mtests <- snp.rhs.tests(cc ~ stratum, family = "binomial", 
     data = subject.support, snp.data = snps.use, tests = blocks)
summary(mtests)


###################################################
### chunk number 37: plot-mtests
###################################################
#line 516 "vignettes/snpStats/inst/doc/snpStats-vignette.Rnw"
pm <- p.value(mtests)
plot(hexbin(-log10(pm), xbin=50))


###################################################
### chunk number 38: qqplot-mtests
###################################################
#line 529 "vignettes/snpStats/inst/doc/snpStats-vignette.Rnw"
qq.chisq(-2 * log(pm), df = 2)