## ----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)