###################################################
### chunk number 1: lkc
###################################################
library(encoDnaseI)
data(rawCD4)
rawCD4


###################################################
### chunk number 2: lkd
###################################################
c19g = rawCD4[ chrnum(19) ]
c19g


###################################################
### chunk number 3: lklk
###################################################
c19gxy = getTrkXY(c19g)
plot(c19gxy)


###################################################
### chunk number 4: doinf
###################################################
clipSnps = function( sms, chrn, lo, hi ) {
 allp = getSnpLocs(sms)
 allp = allp-allp[1] # relative
 ok = allp >= lo & allp <= hi
 thesm = smList(sms)[[1]]
 rsn = colnames(thesm)
 rid = rsn[which(ok)]
 thesm = thesm[, rid, drop=FALSE]
 nn = new.env()
 tmp = list(thesm)
 names(tmp) = as.character(chrn)
 assign("smList", tmp, nn)
 sms@smlEnv = nn
 sms@activeSnpInds=which(ok)
 sms
}

rangeX = function(htrk) {
 range(getTrkXY(htrk)$x)
}


###################################################
### chunk number 5: dogg
###################################################
library(GGtools)
library(GGdata)
if (!exists("hmceuB36")) data(hmceuB36)
rs19g = rangeX(c19g)
h19 = hmceuB36[chrnum(19),]
h19locs = getSnpLocs(hmceuB36[chrnum(19),])[[1]]
goodlocs = which(h19locs[2,] >= rs19g[1] & h19locs[2,] <= rs19g[2])
h19rsn = paste("rs", h19locs[1,goodlocs], sep="")
h19trim = h19[rsid(h19rsn),]
#ok
#c19gf = clipSnps( hmceuB36[chrnum(19),], chrnum(19), rs19g[1], rs19g[2] )
#c19gf


###################################################
### chunk number 6: lkmxi1
###################################################
oo = options()  # don't take warnings on multiple probes... caveat emptor
options(warn=0)
library(GGtools)
showMethods("gwSnpTests")
smxi1 = gwSnpTests(genesym("MXI1")~1-1, h19trim, chrnum(19))
plot(smxi1)
options(oo)


###################################################
### chunk number 7: doj
###################################################
#juxtaPlot = function( trk, ssr ) {
# sy = abs(ssr$trat)
# sx = ssr@locs
# txy = getTrkXY( trk )
# df = data.frame( y = sy, x = sx, type="absT" )
# df = rbind(df, data.frame(y = txy$y, x = txy$x, type = "dnaseI"))
# df$type = factor(df$type, levels=c("dnaseI", "absT"))
# require(lattice)
# xyplot( y~x| type, data=df, layout=c(1,2), 
#    scales=list(relation=list(y="free")))
#}
print(juxtaPlot( c19g, smxi1 ))


###################################################
### chunk number 8: donex
###################################################
oo = options()
options(warn=0)
sOSR2 = gwSnpTests(genesym("OSR2")~1-1, h19trim, chrnum(19))
print(juxtaPlot( c19g, sOSR2 ))
options(oo)


###################################################
### chunk number 9: dow
###################################################
ALICOR( sOSR2, c19g )
ALICOR( smxi1, c19g )


###################################################
### chunk number 10: donsf eval=FALSE
###################################################
## if (interactive()) {
## if (!exists("mads")) mads = apply(exprs(c19gf), 1, mad)
## if (interactive()) fn = featureNames(c19gf)[ which(mads > quantile(mads, .6)) ]
## if (!interactive()) fn = featureNames(c19gf)[ which(mads > quantile(mads, .97)) ]
## n19g = c19gf[ exFeatID(fn), ]
## if (file.exists("tw19g.rda")) load("tw19g.rda")
## if (!exists("tw19g")) tw19g = twSnpScreen( n19g, chr19gmeta, ~., fastAGMfitter)
## if (!file.exists("tw19g.rda")) save(tw19g, file="tw19g.rda")
## if (file.exists("allscor.rda")) load("allscor.rda")
## if (!exists("allscor")) allscor = sapply(tw19g, function(x) {if (inherits(x, "try-error")) return(NA) else return(ALICOR(x, c19g))})
## if (!file.exists("allscor.rda")) save(allscor, file="allscor.rda")
## }
##