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

###################################################
### code chunk number 1: start
###################################################
library(CGEN)


###################################################
### code chunk number 2: load data
###################################################
data(Xdata, package="CGEN")
Xdata[1:5, ]


###################################################
### code chunk number 3: dummy vars
###################################################
for (a in unique(Xdata[, "age.group"])) {
  dummyVar <- paste("age.group_", a, sep="")
  Xdata[, dummyVar] <- 0
  temp <- Xdata[, "age.group"] == a
  if (any(temp)) Xdata[temp, dummyVar] <- 1
}


###################################################
### code chunk number 4: table
###################################################
table(Xdata[, "case.control"], Xdata[, "age.group"], exclude=NULL)


###################################################
### code chunk number 5: snp.logistic
###################################################
mainVars <- c("oral.years", "n.children", "age.group_1",
              "age.group_2", "age.group_3", "age.group_5")
fit <- snp.logistic(Xdata, "case.control", "BRCA.status",
                     main.vars=mainVars, 
                     int.vars=c("oral.years", "n.children"), 
                     strata.var="ethnic.group")


###################################################
### code chunk number 6: summary table
###################################################
getSummary(fit)


###################################################
### code chunk number 7: Wald test
###################################################
getWaldTest(fit, c("BRCA.status", "BRCA.status:oral.years", "BRCA.status:n.children"))


###################################################
### code chunk number 8: snp.list
###################################################
g <- system.file("sampleData", "SNPdata.rda", package="CGEN")
snp.list <- list(file=g, file.type=1, delimiter="|", in.miss="NA")


###################################################
### code chunk number 9: pheno.list
###################################################
f <- system.file("sampleData", "Xdata.txt", package="CGEN")
pheno.list <- list(file=f, id.var="id", file.type=3, delimiter="\t",
              response.var="case.control", strata.var="ethnic.group",
              main.vars=c("oral.years", "n.children"),
              int.vars="oral.years")


###################################################
### code chunk number 10: options list
###################################################
out.file <- "snp.scan.logistic.output.file_h5jb47j.txt"
op <- list(out.file=out.file)


###################################################
### code chunk number 11: snp.scan.logistic
###################################################
ret <- snp.scan.logistic(snp.list, pheno.list, op=op)


###################################################
### code chunk number 12: read data
###################################################
x <- read.table(out.file, sep="\t", header=TRUE)
x[1:5, ]


###################################################
### code chunk number 13: qq plot
###################################################
ret <- QQ.plot(x[, "UML.omnibus.pvalue"])


###################################################
### code chunk number 14: system.file
###################################################
map <- system.file("sampleData", "LocusMapData.txt", package="CGEN")


###################################################
### code chunk number 15: locus map list
###################################################
read.table(map, sep="\t", header=1, nrows=5)
lmap.list <- list(file=map, file.type=3, header=1, delimiter="\t",
                 snp.var="SNP", chrm.var="CHROMOSOME", loc.var="LOCATION")


###################################################
### code chunk number 16: plot setup
###################################################
plot.vars <- c("UML.omnibus.pvalue", "CML.omnibus.pvalue", "EB.omnibus.pvalue")
op <- list(pch=c(21, 22, 23), alt.colors=1)
#plot.new()


###################################################
### code chunk number 17: chromosome.plot
###################################################
ret <- chromosome.plot(out.file, plot.vars, lmap.list, op=op)


###################################################
### code chunk number 18: top hits
###################################################
temp <- sort(x[, "UML.omnibus.pvalue"], index.return=TRUE)$ix
topHits <- x[temp, ][1:10,]
topHits


###################################################
### code chunk number 19: print table
###################################################
table(Xdata$case.control, Xdata$ethnic.group)


###################################################
### code chunk number 20: getMatchedSets
###################################################
library("cluster")
size <- ifelse(Xdata$ethnic.group == 3, 3, 4)
d <- daisy(Xdata[,c("age.group_1","gynSurgery.history","BRCA.history")])
mx <- getMatchedSets(d, CC=TRUE, NN=TRUE, ccs.var = Xdata$case.control, 
	  strata.var = Xdata$ethnic.group, size = size, fixed = TRUE)


###################################################
### code chunk number 21: Xdata
###################################################
mx$CC[1:10]
mx$tblCC
Xdata <- cbind(Xdata, CCStrat = mx$CC, NNStrat = mx$NN)
Xdata[1:5,]


###################################################
### code chunk number 22: snp.matched
###################################################

intVars <- ~ oral.years + n.children
snpVars <- ~ BRCA.status
fit <- snp.matched(Xdata, "case.control", snp.vars=snpVars,
                     main.vars=intVars, int.vars=intVars, 
                     cc.var="CCStrat", nn.var="NNStrat") 


###################################################
### code chunk number 23: summary table 2
###################################################
getSummary(fit, method = c("CLR", "CCL"))


###################################################
### code chunk number 24: Wald test 2
###################################################
getWaldTest(fit$HCL, c( "BRCA.status" , "BRCA.status:oral.years" , "BRCA.status:n.children"))


###################################################
### code chunk number 25: sessionInfo
###################################################
sessionInfo()