snpStats: SnpMatrix and XSnpMatrix classes and methods

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

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


###################################################
### code chunk number 2: select
###################################################
training <- sample(1000, 200)
select <- seq(1, ncol(snps.10),2) 
missing <- snps.10[training, select]
present <- snps.10[training, -select] 
missing 
present 


###################################################
### code chunk number 3: target
###################################################
target <- snps.10[-training, -select]
target


###################################################
### code chunk number 4: imputation-vignette.Rnw:91-93
###################################################
lost <- snps.10[-training, select]
lost


###################################################
### code chunk number 5: positions
###################################################
pos.miss <- snp.support$position[select]
pos.pres <- snp.support$position[-select]


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


###################################################
### code chunk number 7: rule1
###################################################
rules[1:10]


###################################################
### code chunk number 8: rule2
###################################################
rules[c('rs11253563', 'rs2379080')]


###################################################
### code chunk number 9: summary
###################################################
summary(rules)


###################################################
### code chunk number 10: ruleplot
###################################################
plot(rules)


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


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


###################################################
### code chunk number 13: 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")


###################################################
### code chunk number 14: 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")


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


###################################################
### code chunk number 16: imptest-rhs
###################################################
imp2 <- snp.rhs.tests(cc~strata(stratum), family="binomial", 
                      data=subject.support, snp.data=target, 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")


###################################################
### code chunk number 17: impstore
###################################################
imputed <- impute.snps(rules, target, as.numeric=FALSE)


###################################################
### code chunk number 18: uncert1
###################################################
plotUncertainty(imputed[, "rs4880568"])


###################################################
### code chunk number 19: uncert2
###################################################
plotUncertainty(imputed[, "rs2050968"])


###################################################
### code chunk number 20: imptest2
###################################################
imp3 <- single.snp.tests(cc, stratum, data=subject.support,
                        snp.data=imputed, uncertain=TRUE)


###################################################
### code chunk number 21: imp3
###################################################
imp3[1:5]
imp[1:5]


###################################################
### code chunk number 22: mach
###################################################
path <- system.file("extdata/mach1.out.mlprob.gz", package="snpStats")
mach <- read.mach(path)
plotUncertainty(mach[,50])


###################################################
### code chunk number 23: class-imp-obs
###################################################
class(imp)


###################################################
### code chunk number 24: 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=target, rules=rules, score=TRUE)


###################################################
### code chunk number 25: imputation-vignette.Rnw:345-347
###################################################
class(obs)
class(imp)


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


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


###################################################
### code chunk number 28: sign
###################################################
table(effect.sign(obs))


###################################################
### code chunk number 29: switch
###################################################
effect.sign(obs)[1:6]
sw.obs <- switch.alleles(obs, 1:3)
class(sw.obs)
effect.sign(sw.obs)[1:6]