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


###################################################
### chunk number 2: select
###################################################
sel <- seq(1, ncol(snps.10),2) 
missing <- snps.10[,sel]
present <- snps.10[,-sel] 
missing 
present 


###################################################
### chunk number 3: positions
###################################################
pos.miss <- snp.support$position[sel]
pos.pres <- snp.support$position[-sel]


###################################################
### chunk number 4: rules
###################################################
rules <- snp.imputation(present, missing, pos.pres, pos.miss)


###################################################
### chunk number 5: rule1
###################################################
rules[1:10]


###################################################
### chunk number 6: summary
###################################################
summary(rules)


###################################################
### chunk number 7: ruleplot
###################################################
plot(rules)


###################################################
### chunk number 8: imptest
###################################################
imp <- single.snp.tests(cc, stratum, data=subject.support,
                        snp.data=present, rules=rules)


###################################################
### chunk number 9: realtest
###################################################
obs <- single.snp.tests(cc, stratum, data=subject.support, snp.data=missing)


###################################################
### chunk number 10: compare
###################################################
logP.imp <- -log10(p.value(imp, df=1))
logP.obs <- -log10(p.value(obs, df=1))
hb <- hexbin(logP.obs, logP.imp, xbin=50)
sp <- plot(hb)
hexVP.abline(sp$plot.vp, 0, 1, col="black")


###################################################
### chunk number 11: best
###################################################
use <- imputation.r2(rules)>0.9
hb <- hexbin(logP.obs[use], logP.imp[use], xbin=50)
sp <- plot(hb)
hexVP.abline(sp$plot.vp, 0, 1, col="black")


###################################################
### chunk number 12: rsqmaf
###################################################
hb <- hexbin(imputation.maf(rules), imputation.r2(rules), xbin=50)
sp <- plot(hb)


###################################################
### chunk number 13: imptest-rhs
###################################################
imp2 <- snp.rhs.tests(cc~strata(stratum), family="binomial", 
                      data=subject.support, snp.data=present, rules=rules)
logP.imp2 <- -log10(p.value(imp2))
hb <- hexbin(logP.obs, logP.imp2, xbin=50)
sp <- plot(hb)
hexVP.abline(sp$plot.vp, 0, 1, col="black")


###################################################
### chunk number 14: class-imp-obs
###################################################
class(imp)
class(obs)


###################################################
### chunk number 15: save-scores
###################################################
obs <- single.snp.tests(cc, stratum, data=subject.support, snp.data=missing, 
                        score=TRUE)
imp <- single.snp.tests(cc, stratum, data=subject.support,
                        snp.data=present, rules=rules, score=TRUE)
class(obs)
class(imp)


###################################################
### chunk number 16: pool
###################################################
both <- pool(obs, imp)
class(both)
both[1:5]


###################################################
### chunk number 17: pool-score
###################################################
both <- pool(obs, imp, score=TRUE)
class(both)


###################################################
### chunk number 18: sign
###################################################
table(effect.sign(obs))


###################################################
### chunk number 19: switch
###################################################
effect.sign(obs)[1:6]
sw.obs <- switch.alleles(obs, c("rs12773042", "rs10904596"))
class(sw.obs)
effect.sign(sw.obs)[1:6]