###################################################
### chunk number 1: customizeSweave
###################################################
#line 59 "/home/schillert/talks/2010_12-Floripa-IBC-Andreas/tutorial/scripts/final/illustration-GWA_data_handling.Rnw"
options(width = 65, prompt = "R> ", continue = "   ")


###################################################
### chunk number 2: loadData
###################################################
#line 85 "/home/schillert/talks/2010_12-Floripa-IBC-Andreas/tutorial/scripts/final/illustration-GWA_data_handling.Rnw"
library(GenABEL)                     

pheFile <- 'wgDat.phe'
rawFile <- 'wgDat.raw'

wgDat <- load.gwaa.data(pheFile, rawFile)


###################################################
### chunk number 3: 
###################################################
#line 99 "/home/schillert/talks/2010_12-Floripa-IBC-Andreas/tutorial/scripts/final/illustration-GWA_data_handling.Rnw"
strOptions(strict.width = "yes")


###################################################
### chunk number 4: strGwaa
###################################################
#line 103 "/home/schillert/talks/2010_12-Floripa-IBC-Andreas/tutorial/scripts/final/illustration-GWA_data_handling.Rnw"
str(wgDat)


###################################################
### chunk number 5: tabAff
###################################################
#line 108 "/home/schillert/talks/2010_12-Floripa-IBC-Andreas/tutorial/scripts/final/illustration-GWA_data_handling.Rnw"
table(wgDat@phdata$affection)


###################################################
### chunk number 6: recodeAff
###################################################
#line 115 "/home/schillert/talks/2010_12-Floripa-IBC-Andreas/tutorial/scripts/final/illustration-GWA_data_handling.Rnw"
wgDat@phdata$aff.01 <- wgDat@phdata$affection - 1
with(wgDat@phdata, table(affection, aff.01))


###################################################
### chunk number 7: 
###################################################
#line 140 "/home/schillert/talks/2010_12-Floripa-IBC-Andreas/tutorial/scripts/final/illustration-GWA_data_handling.Rnw"
idSummar <- perid.summary(wgDat)
head(idSummar)


###################################################
### chunk number 8: 
###################################################
#line 150 "/home/schillert/talks/2010_12-Floripa-IBC-Andreas/tutorial/scripts/final/illustration-GWA_data_handling.Rnw"
hetMean <- mean(idSummar$Het)
hetSd <- sd(idSummar$Het)
hetThreshLow <- hetMean - 3 * hetSd
hetThreshUpp <- hetMean + 3 * hetSd

removeIdx <- with(idSummar, which(CallPP < 0.97 | Het < hetThreshLow | Het > hetThreshUpp))
idSummar$keep <- TRUE
idSummar$keep[removeIdx] <- FALSE
keepIDs <- row.names(idSummar[idSummar$keep, ])

wgDatIdClean <- wgDat[keepIDs, ]


###################################################
### chunk number 9: call_vs_het
###################################################
#line 169 "/home/schillert/talks/2010_12-Floripa-IBC-Andreas/tutorial/scripts/final/illustration-GWA_data_handling.Rnw"
with(idSummar,  plot(CallPP, Het, pch = 19, ##xlim = c(0.95, 1),
                     xlab = "Call fraction", ylab = "Proportion of heterozygosity",
                     col = ifelse(keep, "black", "red")))
abline(h = hetThreshLow, col = "lightgrey")
abline(h = hetThreshUpp, col = "lightgrey")
abline(v = 0.97, col = "lightgrey")


###################################################
### chunk number 10: pSim
###################################################
#line 201 "/home/schillert/talks/2010_12-Floripa-IBC-Andreas/tutorial/scripts/final/illustration-GWA_data_handling.Rnw"
pwSim <- ibs(wgDatIdClean)
pwDist <- as.dist(0.5 - pwSim)
mdsDat <- cmdscale(pwDist)
plot(mdsDat[, 1], mdsDat[, 2], xlab = "Component 1",
     ylab = "Component 2", pch = 19)


###################################################
### chunk number 11: snpSum
###################################################
#line 234 "/home/schillert/talks/2010_12-Floripa-IBC-Andreas/tutorial/scripts/final/illustration-GWA_data_handling.Rnw"
casesIDs <- subset(wgDatIdClean@phdata, affection == 2, id, drop = TRUE)
controlsIDs <- subset(wgDatIdClean@phdata, affection == 1, id, drop = TRUE)
sumMaf <- summary(wgDat)$Q.2
maf <- data.frame(maf = pmin(sumMaf, 1- sumMaf))

sumControls <- summary(wgDatIdClean[controlsIDs, ])[, c("Pexact", "CallRate")]
names(sumControls) <- c("pHWE", "cfControls")
sumCases <- summary(wgDatIdClean[casesIDs, ])[, "CallRate", drop = FALSE]
names(sumCases) <- "cfCases"
snpSummar <- do.call(cbind, list(maf, sumControls, sumCases))

snpRemoveIdx <- with(snpSummar, which( maf < 0.01 | pHWE < 0.0001 | cfCases < 0.98 | cfControls < 0.98))

snpSummar$keep <- TRUE
snpSummar$keep[snpRemoveIdx] <- FALSE
keepSNPs <- row.names(snpSummar[snpSummar$keep, ])

wgDatClean <- wgDatIdClean[, keepSNPs]


###################################################
### chunk number 12: venn
###################################################
#line 258 "/home/schillert/talks/2010_12-Floripa-IBC-Andreas/tutorial/scripts/final/illustration-GWA_data_handling.Rnw"
source("plot-VennDiagram.R")


###################################################
### chunk number 13: venn
###################################################
#line 267 "/home/schillert/talks/2010_12-Floripa-IBC-Andreas/tutorial/scripts/final/illustration-GWA_data_handling.Rnw"
mafSNPs <- row.names(subset(snpSummar, maf < 0.01))
hweSNPs <- row.names(subset(snpSummar, pHWE < 0.0001))
crSNPs <- row.names(subset(snpSummar, cfCases < 0.98 | cfControls <0.98))

qcVenn(mafSNPs, hweSNPs, crSNPs, labels = c("MAF", "HWE", "CF"),
       numberSnps = nrow(snpSummar))


###################################################
### chunk number 14: mlreg
###################################################
#line 294 "/home/schillert/talks/2010_12-Floripa-IBC-Andreas/tutorial/scripts/final/illustration-GWA_data_handling.Rnw"
assocRes <-  mlreg(aff.01 ~ 1, data = wgDatClean, 
                   gtmode = "additive", trait.type  = "binomial") 
plot(assocRes, main = "", ystart = 2) 


###################################################
### chunk number 15:  eval=FALSE
###################################################
## #line 310 "/home/schillert/talks/2010_12-Floripa-IBC-Andreas/tutorial/scripts/final/illustration-GWA_data_handling.Rnw"
## help("scan.gwaa-class")


###################################################
### chunk number 16: 
###################################################
#line 319 "/home/schillert/talks/2010_12-Floripa-IBC-Andreas/tutorial/scripts/final/illustration-GWA_data_handling.Rnw"
toLatex(sessionInfo(), locale = FALSE)