GEEs

## ----style, echo = FALSE, results = 'asis'------------------------------------
BiocStyle::markdown()

## ----dodisp-------------------------------------------------------------------
library(MASS)
library(BatchJobs)
data(anorexia)  # N = 72
myr = makeRegistry("abc", file.dir=tempfile())
chs = chunk(1:nrow(anorexia), n.chunks=18) # 4 recs/chunk
f = function(x) anorexia[x,]
options(BBmisc.ProgressBar.style="off")
batchMap(myr, f, chs)
showStatus(myr)

## ----execdisp, results="hide", echo=FALSE-------------------------------------
suppressMessages(submitJobs(myr,progressbar=FALSE))
waitForJobs(myr) 

## ----fakk, eval=FALSE---------------------------------------------------------
#  submitJobs(myr)
#  waitForJobs(myr)

## ----lkres--------------------------------------------------------------------
loadResult(myr,1)

## ----dopar--------------------------------------------------------------------
library(parglms)
pp = parGLM( Postwt ~ Treat + Prewt, myr,
   family=gaussian, binit = c(0,0,0,0), maxit=10, tol=.001 )
names(pp)
pp$coef

## ----defineUpdate, eval=FALSE-------------------------------------------------
#  litdec = function(grWithFDR) {
#   tmp = grWithFDR
#   library(gQTLstats)
#   if (!exists("hmm878")) data(hmm878)
#   seqlevelsStyle(hmm878) = "NCBI"
#   library(GenomicRanges)
#   ov = findOverlaps(tmp, hmm878)
#   states = hmm878$name
#   states[ which(states %in% c("13_Heterochrom/lo", "14_Repetitive/CNV",
#        "15_Repetitive/CNV")) ] = "hetrep_1315"
#   levs = unique(states)
#   tmp$hmmState = factor(rep("hetrep_1315", length(tmp)),levels=levs)
#   tmp$hmmState = relevel(tmp$hmmState, "hetrep_1315")
#   tmp$hmmState[ queryHits(ov) ] = factor(states[ subjectHits(ov) ],
#     levels=levs)
#   if (!exists("gwrngs19")) data(gwrngs19)
#   library(GenomeInfoDb)
#   seqlevelsStyle(gwrngs19) = "NCBI"
#   tmp$isGwasHit = 1*(tmp %in% gwrngs19)
#   tmp
#  }
#  
#  decorate = function(grWithFDR) {
#  #
#  # the data need a distance/MAF filter
#  #
#   library(gQTLstats)
#   data(filtFDR)
#   if (!exists("hmm878")) data(hmm878)
#   library(gwascat)
#   if (!exists("gwrngs19")) data(gwrngs19)
#   if (!exists("gwastagger")) data(gwastagger) # will use locations here
#   library(GenomeInfoDb)
#   seqlevelsStyle(hmm878) = "NCBI"
#   seqlevelsStyle(gwrngs19) = "NCBI"
#   seqlevelsStyle(gwastagger) = "NCBI"
#   tmp = grWithFDR
#   tmp$isGwasHit = 1*(tmp %in% gwrngs19)
#   tmp$isGwasTagger = 1*(tmp %in% gwastagger)
#   #levs = unique(hmm878$name)
#   library(GenomicRanges)
#   ov = findOverlaps(tmp, hmm878)
#   states = hmm878$name
#   states[ which(states %in% c("13_Heterochrom/lo", "14_Repetitive/CNV",
#        "15_Repetitive/CNV")) ] = "hetrep_1315"
#   levs = unique(states)
#   tmp$hmmState = factor(rep("hetrep_1315", length(tmp)),levels=levs)
#   tmp$hmmState = relevel(tmp$hmmState, "hetrep_1315")
#   tmp$hmmState[ queryHits(ov) ] = factor(states[ subjectHits(ov) ],
#     levels=levs)
#   tmp$estFDR = getFDRfunc(filtFDR)( tmp$chisq )
#   tmp$fdrcat = cut(tmp$estFDR, c(-.01, .01, .05, .1, .25, .5, 1.01))
#   tmp$fdrcat = relevel(tmp$fdrcat, "(0.5,1.01]")
#   #tmp$distcat = cut(tmp$mindist, c(-1,0,1000,5000,10000,50000,100000,250000,500001))
#   tmp$distcat = cut(tmp$mindist, c(-1,0,1000,5000,10000,25000,50001))
#   #tmp$distcat = relevel(tmp$distcat, "(2.5e+05,5e+05]")
#   tmp$distcat = relevel(tmp$distcat, "(2.5e+04,5e+04]")
#   tmp$MAFcat = cut(tmp$MAF, c(.049, .075, .1, .25, .51))
#   tmp$MAFcat = relevel(tmp$MAFcat, "(0.25,0.51]")
#   kp = c("seqnames", "start", "probeid", "snp", "estFDR", "fdrcat", "hmmState",
#    "distcat", "MAFcat", "isGwasHit", "isGwasTagger")
#   names(tmp) = NULL
#   as(tmp, "data.frame")[,kp]
#  }

## ----doge, eval=FALSE---------------------------------------------------------
#  suppressPackageStartupMessages({
#  library(geuvStore2)
#  library(gQTLBase)
#  library(gQTLstats)
#  })
#  prst = g17transRegistry()

## ----domod,eval=FALSE---------------------------------------------------------
#  prst$extractor = function(store,i) litdec(loadResult(store,i)[[1]])
#  p1 = parGLM( isGwasHit ~ hmmState, prst,
#     family=binomial, binit=rep(0,13), tol=.001,
#     maxit = 10 )
#  summaryPG(p1)
#  #ans= list(coef=p1$coef, s.e.=sqrt(diag(p1$eff.var)))
#  #ans$z = ans[[1]]/ans[[2]]
#  #do.call(cbind, ans)