### R code from vignette source 'vignettes/GGtools/inst/doc/GGtools_2011.Rnw'

###################################################
### code chunk number 1: lkd
###################################################
library(GGtools)
if (!exists("hmceuB36.2021")) hmceuB36.2021 <- getSS("GGtools", c("20", "21"))
hmceuB36.2021


###################################################
### code chunk number 2: lkn
###################################################
library("illuminaHumanv1.db")
cpn = get("CPNE1", revmap(illuminaHumanv1SYMBOL))
use = intersect(cpn, featureNames(hmceuB36.2021))
if (length(use) == 0) stop("probe not on array")
use = use[1]


###################################################
### code chunk number 3: trim
###################################################
hmlit = hmceuB36.2021[probeId(use),]
hmlit = hmlit[ chrnum("20"), ]
hmlit


###################################################
### code chunk number 4: doe
###################################################
tname = function() gsub(".*/|.*\\\\", "", tempfile())
e1 = eqtlTests(hmlit, ~male, targdir=tname(), runname=tname())
e1
class(e1)


###################################################
### code chunk number 5: doq
###################################################
e1[rsid("rs6060535"), probeId("GI_23397697-A")]


###################################################
### code chunk number 6: dotop
###################################################
topFeats(probeId(cpn), mgr=e1, ffind=1)


###################################################
### code chunk number 7: get09
###################################################
library(SNPlocs.Hsapiens.dbSNP.20090506)
c20l = getSNPlocs("chr20")
c20 = GRanges(IRanges(c20l$loc, width=1), seqnames="chr20")
rs20 = paste("rs", c20l$RefSNP_id, sep="")
names(c20) = rs20
length(c20)
c20[1:3,]


###################################################
### code chunk number 8: lkn
###################################################
tn = rownames(e1@fflist[[1]])
tn[1:10]


###################################################
### code chunk number 9: lkman
###################################################
manhPlot(cpn, mgr=e1, ffind=1, locGRanges=c20, cex=3)


###################################################
### code chunk number 10: dook
###################################################
okn = intersect(tn, names(c20))
c20ok = c20[okn]


###################################################
### code chunk number 11: getp
###################################################
chis1 = e1[, probeId(cpn)][[1]][okn,]
ml10p = -log10(1-pchisq(chis1, 1))


###################################################
### code chunk number 12: dob
###################################################
elementMetadata(c20ok)$score = ml10p
c20ok[1:5,]


###################################################
### code chunk number 13: doex
###################################################
library(rtracklayer)
export(c20ok, "cpne1.wig")


###################################################
### code chunk number 14: lkgsc
###################################################
top200 = c20ok[order(elementMetadata(c20ok)$score, decreasing=TRUE)[1:200]]
top200[196:200]


###################################################
### code chunk number 15: gettab (eval = FALSE)
###################################################
## library(rtracklayer)
## ss = browserSession()
## genome(ss) = "hg18"
## GG = GRanges(seqnames="chr20", IRanges(30000000,40000000))
## oq = ucscTableQuery(ss, "oreganno", range=GG)
## toq = track(oq)
## ov = subsetByOverlaps(toq, top200)
## ov


###################################################
### code chunk number 16: domore (eval = FALSE)
###################################################
## tableName(oq) = "oregannoAttr"
## featdat = getTable(oq)
## featdat[ featdat$id %in% ov$name, ]


###################################################
### code chunk number 17: lke
###################################################
NG = length(featureNames(hmceuB36.2021))
hmct = hmceuB36.2021[ sample(1:NG, size=20, replace=FALSE), ]
applier = lapply
chkmult = function() {length(grep("^multicore$", installed.packages()[,1])) > 0}
if (chkmult()) {
  library(multicore)
  options(cores=min(c(multicore:::volatile$detectedCores, 6)))
  applier = mclapply
  }
unix.time(e2 <- eqtlTests(hmct, ~male, geneApply=applier, targdir=tempdir(), runname=tname()))


###################################################
### code chunk number 18: lkff
###################################################
e2@fflist[[1]]
sapply(e2@fflist,dim)


###################################################
### code chunk number 19: doext (eval = FALSE)
###################################################
## edname = paste("exdem", tname(), sep="")
## if (interactive() & file.exists(edname)) stop(paste("will not remove existing '", edname, "' folder"))
## externalize(hmceuB36.2021, edname)
## tar(tna <- paste(edname,".tar.gz", sep=""), edname, compression="gzip")
## install.packages(tna, repos=NULL, type="source")
## try(system(paste("rm -rf", tna)))
## try(system(paste("rm -rf", edname)))


###################################################
### code chunk number 20: lkd (eval = FALSE)
###################################################
## gc()
## s1 = getSS(edname, "20")
## gc()
## rm(s1)
## gc()


###################################################
### code chunk number 21: mymap
###################################################
mymap = lapply(c("20", "21"), function(x) get(x, revmap(illuminaHumanv1CHR))[1:20])
mymap = lapply(mymap, function(x) intersect(x, featureNames(hmceuB36.2021)))
names(mymap) = c("20", "21")
unix.time(md <-  makeDiagDirector("GGtools", mymap, rhs=~male, sleeplen=2, runname=tname(),
    targdir=tname()))
md


###################################################
### code chunk number 22: getb
###################################################
ps = lapply(md@mgrs, function(x) colnames(x@fflist[[1]]))
names(ps) = names(md@mgrs)
ps


###################################################
### code chunk number 23: getgr
###################################################
getgr = function(x,map) abs(sapply(mget(x, map), "[", 1))
st = lapply(ps, getgr, illuminaHumanv1CHRLOC)
en = lapply(ps, getgr, illuminaHumanv1CHRLOCEND)
sp = names(ps)
grl = list()
for (i in 1:length(ps)) {
  grl[[i]] = GRanges(seqnames=sp[i], IRanges(st[[i]], en[[i]]))
  names(grl[[i]]) = ps[[i]]
  }
names(grl) = names(ps)


###################################################
### code chunk number 24: getsr
###################################################
seqlevels(c20) <- sub("chr", "", seqlevels(c20))
c21l = getSNPlocs("chr21")  # hg18
c21 = GRanges(IRanges(c21l$loc, width=1), seqnames="21")
rs21 = paste("rs", c21l$RefSNP_id, sep="")
names(c21) = rs21
srl = list(`20`=c20, `21`=c21)


###################################################
### code chunk number 25: lksc
###################################################
csc = cisProxScores( direc=md, dradset = c(0,50000,100000, 500000), snpGRL = srl, geneGRL = grl,
  geneApply=lapply, ffind=1 )


###################################################
### code chunk number 26: lklap
###################################################
lapply(csc, lapply, sapply,  function(x)max(unlist(x)))


###################################################
### code chunk number 27: lks
###################################################
sessionInfo()