### R code from vignette source 'tut11.Rnw'

###################################################
### code chunk number 1: dogg
###################################################
intsave = function(...) NULL # could set as save to speed runs
date()
library(GGtools)
library(GGdata)
c17 = getSS("GGdata", "17", renameChrs="chr17")
class(c17)  # smlSet links SnpMatrix instances and expression data
c17


###################################################
### code chunk number 2: dote
###################################################
t1 = gwSnpTests(genesym("CHRNE")~male, c17, chrnum("chr17"))
topSnps(t1)


###################################################
### code chunk number 3: lkass
###################################################
par(mfrow=c(2,2))
plot_EvG(genesym("CHRNE"), rsid("rs16954243"), c17)
plot_EvG(genesym("CHRNE"), rsid("rs7214776"), c17)
plot_EvG(genesym("CHRNE"), rsid("rs2302321"), c17)
plot_EvG(genesym("CHRNE"), rsid("rs8070572"), c17)
par(mfrow=c(1,1))


###################################################
### code chunk number 4: getresp
###################################################
chrneExpr = as.numeric(exprs(c17[genesym("CHRNE"),]))
summary(chrneExpr)


###################################################
### code chunk number 5: getg
###################################################
num243 = as( smList(c17)[["chr17"]][, "rs16954243"], "numeric" )
table(num243)


###################################################
### code chunk number 6: lklm
###################################################
summary(lm(chrneExpr~num243+male, data=pData(c17)))


###################################################
### code chunk number 7: lkactual
###################################################
snp.rhs.tests(chrneExpr~1, data=pData(c17), 
   snp.data=smList(c17)[["chr17"]][, "rs16954243"],
   family="gaussian")


###################################################
### code chunk number 8: lkva
###################################################
var.test(chrneExpr[num243==0], chrneExpr[num243==1])


###################################################
### code chunk number 9: lktu
###################################################
library(ggtut)
f1 = observed17ceu()
f1
f1@call


###################################################
### code chunk number 10: lkpr
###################################################
topFeats(probeId("GI_44662825-I"), mgr=f1, ffind=1)


###################################################
### code chunk number 11: dosurv
###################################################
options(digits=4)
bestApply = lapply
if ("multicore" %in% installed.packages()[,1]) {
  library(multicore)
  options(cores=10)  # check!
  bestApply = mclapply
}
allpro = probesManaged(f1,1)
if (!exists("tops")) tops = bestApply( allpro, 
  function(x) topFeats(probeId(x), mgr=f1, ffind=1, n=5 ) )
names(tops) = allpro


###################################################
### code chunk number 12: lkp
###################################################
lapply(tops[1:4], function(x)1-pchisq(x,1))


###################################################
### code chunk number 13: lkperm
###################################################
permf1 = onePerm17ceu()
if (!exists("permtops")) permtops = bestApply( allpro, 
   function(x) topFeats(probeId(x),
   mgr=permf1, ffind=1, n=5))
names(permtops) = allpro
lapply(permtops[1:4], function(x)1-pchisq(x,1))


###################################################
### code chunk number 14: maxass
###################################################
maxassoc = sapply(tops, "[", 1)
maxassoc[1:5]


###################################################
### code chunk number 15: maxpass
###################################################
maxaperm = sapply(permtops, "[", 1)
maxaperm[1:5]


###################################################
### code chunk number 16: lk99
###################################################
p99 = quantile(maxaperm, .99)
p99
sum(maxassoc > p99)


###################################################
### code chunk number 17: getg
###################################################
haseqtl = which(maxassoc > p99)
pweq = allpro[haseqtl]
unlist(mget(pweq, illuminaHumanv1SYMBOL))


###################################################
### code chunk number 18: getsn
###################################################
tmp = names(maxassoc)[haseqtl]
tmp = gsub(".rs", "%.rs", tmp)
pids = sapply(strsplit(tmp, "%."), "[", 1)
rsids = sapply(strsplit(tmp, "%."), "[", 2)
par(mfrow=c(2,2))
for (i in 1:4) plot_EvG( probeId(pids[i]),
   rsid(rsids[i]), c17)
par(mfrow=c(1,1))


###################################################
### code chunk number 19: getlocs
###################################################
data(snpgr17)
length(snpgr17)
snpgr17[1:3]
length(intersect(names(snpgr17), snpsManaged(f1,1)))


###################################################
### code chunk number 20: lkdb
###################################################
library(GenomicFeatures)
if (!exists("txdb")) 
    txdb = loadFeatures(system.file("sqlite/hg18.txdb.sqlite", 
      package="ggtut"))
txg = transcriptsBy(txdb, "gene")
txg[1:3]


###################################################
### code chunk number 21: mkg17 (eval = FALSE)
###################################################
## pn = probesManaged(f1,1)
## library(illuminaHumanv1.db)
## pn.eg = unlist(mget(pn, illuminaHumanv1ENTREZID))
## pn.eg = na.omit(pn.eg)
## eg.pn = names(pn.eg)
## names(eg.pn) = pn.eg
## txg17 = txg[ intersect(names(txg), pn.eg) ]
## extents = function(x) {
##   y = x[seqnames(x)=="chr17"]; c(min(start(y)),max(end(y)))
##   }  # watch for random
## ssnr = lapply( txg17, function(z) try(extents(z)) )
## firsts = sapply(ssnr, function(x) {if(is.finite(x[1])) return(x[1]); NA})
## if (any(is.na(firsts))) ssnr = ssnr[-which(is.na(firsts))]
## firsts = sapply(ssnr, function(x) {if(is.numeric(x[1])) return(x[1]); NA})
## lasts = sapply(ssnr, function(x) {if(is.numeric(x[2])) return(x[2]); NA})
## g17rngsnr = GRanges(seqnames="chr17", 
##     IRanges(firsts,lasts), probeid=eg.pn[names(ssnr)])


###################################################
### code chunk number 22: dodat
###################################################
data(g17rngsnr)
g17rngsnr


###################################################
### code chunk number 23: docps (eval = FALSE)
###################################################
## df1 = new("multiCisDirector", mgrs=list(obs17=f1))
## if (!exists("CPS17")) {
## CPS17 = cisProxScores(dradset=c(50000,2e6), direc=df1, 
##     snpGRL=list(obs17=snpgr17), geneGRL=list(obs17=g17rngsnr), ffind=1)
##     }
## permdf1 = new("multiCisDirector", mgrs=list(obs17=permf1))
## if (!exists("PERMCPS17")) {
## PERMCPS17 = cisProxScores(dradset=c(50000,2e6), direc=permdf1, 
##     snpGRL=list(obs17=snpgr17), geneGRL=list(obs17=g17rngsnr), ffind=1)
##     }


###################################################
### code chunk number 24: dodatcps
###################################################
if (!exists("CPS17")) data(CPS17)
if (!exists("PERMCPS17")) data(PERMCPS17)


###################################################
### code chunk number 25: lksb
###################################################
sb1 = scoresByGenes(CPS17, as.GRanges=FALSE)
lapply(sb1[1:3], "[", 1:10)
summary(sapply(sb1,length))
permsb1 = scoresByGenes(PERMCPS17, as.GRanges=FALSE)
summary(unlist(sb1))
summary(unlist(permsb1))


###################################################
### code chunk number 26: dos
###################################################
rsb1 = scoresByGenes(CPS17, as.GRanges=TRUE, snpGR=snpgr17,
    scoreConverter=function(x)-log10(1-pchisq(x,1)))
prsb1 = scoresByGenes(PERMCPS17, as.GRanges=TRUE, snpGR=snpgr17,
    scoreConverter=function(x)-log10(1-pchisq(x,1)))
library(rtracklayer)
export(rsb1, "obs17.wig")
export(prsb1, "perm17.wig")


###################################################
### code chunk number 27: dobigger (eval = FALSE)
###################################################
## rsb2 = scoresByGenes(CPS17, intvind=2, as.GRanges=TRUE, snpGR=snpgr17,
##     scoreConverter=function(x)-log10(1-pchisq(x,1)))
## export(rsb2, "obs17p2m.wig")


###################################################
### code chunk number 28: ggg
###################################################
p99bysnp = quantile(elementMetadata(prsb1)$score, .99)
sum(elementMetadata(rsb1)$score > p99bysnp)
rsb1hi = rsb1[ which(elementMetadata(rsb1)$score > p99bysnp)]
library(ChIPpeakAnno)
data(TSS.human.NCBI36)
TSS.human.NCBI36[1:3,]
rsb1hi = as(rsb1hi, "RangedData")
rsb1report = annotatePeakInBatch(rsb1hi, AnnotationData=TSS.human.NCBI36)
rsb1report[1:4,]
table(rsb1report$insideFeature)


###################################################
### code chunk number 29: dodistpl
###################################################
plot(rsb1hi$score ~ rsb1report$distancetoFeature)


###################################################
### code chunk number 30: getr
###################################################
data(rules.n43)
rules.n43[1:5]
summary(rules.n43)


###################################################
### code chunk number 31: dom (eval = FALSE)
###################################################
## tempfolder = function () 
## {
##     z = tempfile()
##     system(paste("mkdir", z))
##     z
## }
## obsfold = tempfolder()
## p2keep = probesManaged(f1,1)
## c17 = getSS("GGdata", "17", renameChrs="chr17", wrapperEndo=dropMonomorphies,
##    probesToKeep=p2keep)
## if (!exists("rf1")) {
##  if (file.exists("rf1.rda")) load("rf1.rda") else {
##     rf1 = ieqtlTests( c17, ~male, targdir=obsfold,
##       geneApply=bestApply, shortfac=10,
##       rules=rules.n43 )
##     intsave(rf1, file="rf1.rda")
##     }
## }


###################################################
### code chunk number 32: lkmmrevm (eval = FALSE)
###################################################
## newsn = snpsManaged(rf1,1)
## extSNP = newsn[grep("chr17:", newsn)]
## elocs = as.numeric(gsub("chr17:", "", extSNP))
## newr = GRanges(seqnames="chr17", IRanges(elocs,width=1))
## names(newr) = extSNP
## extsnpgr17 = c(snpgr17, newr)
## rdf1 = new("multiCisDirector", mgrs=list(imp17=rf1))
## if (file.exists("IMP17.rda")) load("IMP17.rda")
## if (!exists("IMP17")) {
## IMP17 = cisProxScores(dradset=c(50000,2e6), direc=rdf1,
##     snpGRL=list(imp17=extsnpgr17), geneGRL=list(imp17=g17rngsnr), ffind=1)
## intsave(IMP17, file="IMP17.rda")
##     }
## rsb2 = scoresByGenes(IMP17, as.GRanges=TRUE, snpGR=extsnpgr17,
##     scoreConverter=function(x)-log10(1-pchisq(x,1)))
## rsb3 = scoresByGenes(IMP17, as.GRanges=TRUE, snpGR=extsnpgr17, intvind=2,
##     scoreConverter=function(x)-log10(1-pchisq(x,1)))


###################################################
### code chunk number 33: doimpperm (eval = FALSE)
###################################################
## perfold = tempfolder()
## if (!exists("rf1_perm")) {
##  if (file.exists("rf1_perm.rda")) load("rf1_perm.rda") else {
##     rf1_perm = ieqtlTests( permEx(c17), ~male, targdir=perfold,
##       geneApply=bestApply, shortfac=10,
##       rules=rules.n43 )
##     intsave(rf1_perm, file="rf1_perm.rda")
##     }
## }
## rdf1_perm = new("multiCisDirector", mgrs=list(imp17=rf1_perm))
## if (file.exists("IMP17_PERM.rda")) load("IMP17_PERM.rda")
## if (!exists("IMP17_PERM")) {
## IMP17_PERM = cisProxScores(dradset=c(50000,2e6), direc=rdf1_perm,
##     snpGRL=list(imp17=extsnpgr17), geneGRL=list(imp17=g17rngsnr), ffind=1)
## intsave(IMP17_PERM, file="IMP17_PERM.rda")
## }
## rsb2_perm = scoresByGenes(IMP17_PERM, as.GRanges=TRUE, snpGR=extsnpgr17,
##     scoreConverter=function(x)-log10(1-pchisq(x,1)))
## rsb2_list = scoresByGenes(IMP17, as.GRanges=FALSE)
## rsb2_perm_list = scoresByGenes(IMP17_PERM, as.GRanges=FALSE)


###################################################
### code chunk number 34: doq
###################################################
data(rsb2_perm_list)
p99i = quantile( sapply(rsb2_perm_list, max), .99 )
data(rsb2_list)
sum( sapply( rsb2_list, max) > p99i )


###################################################
### code chunk number 35: doq
###################################################
data(permsb1)
p99u = quantile( sapply(permsb1, max), .99 )
data(sb1)
sum( sapply( sb1, max) > p99u )


###################################################
### code chunk number 36: lkt
###################################################
length(unlist(rsb2_list))


###################################################
### code chunk number 37: lktt
###################################################
length(unlist(sb1))


###################################################
### code chunk number 38: srchcd
###################################################
library(GGtools)
c17 = getSS("GGdata", "17")
get("CD79B", revmap(illuminaHumanv1SYMBOL))


###################################################
### code chunk number 39: domocd
###################################################
lkcd1 = gwSnpTests(probeId("GI_11038675-A")~male, c17, chrnum("17")) 
topSnps(lkcd1)


###################################################
### code chunk number 40: getpr
###################################################
pct = prcomp(t(exprs(c17)))


###################################################
### code chunk number 41: lkpc
###################################################
plot(pct)


###################################################
### code chunk number 42: addsm
###################################################
DF = data.frame(pct$x[,1:4])
pData(c17) = cbind(pData(c17), DF)


###################################################
### code chunk number 43: docd2
###################################################
lkcd2 = gwSnpTests(probeId("GI_11038675-A")~PC1+PC2+PC3+PC4+male, c17, chrnum("17")) 
topSnps(lkcd2)


###################################################
### code chunk number 44: getdrops
###################################################
table(nsn <- as(smList(c17)[[1]][,"rs2584597"], "numeric"))
drop = which(is.na(nsn))


###################################################
### code chunk number 45: newstuff (eval = FALSE)
###################################################
## mod = model.matrix(~as(smList(c17)[[1]][,"rs2584597"], "numeric"))
## mod0 = model.matrix(~1, data=pData(c17[,-drop]))
## library(sva)
## if (!exists("SVA1") & file.exists("SVA1.rda")) load("SVA1.rda") else {
##    SVA1 = sva(exprs(c17[,-drop]), mod, mod0); 
##    intsave(SVA1, file="SVA1.rda"); 
## }


###################################################
### code chunk number 46: getsvas
###################################################
if (!exists("SVA1")) data(SVA1)


###################################################
### code chunk number 47: lknv
###################################################
SVA1$n.sv


###################################################
### code chunk number 48: domote
###################################################
SVDF = data.frame(SVA1$sv)
c17d = c17[,-drop]
pData(c17d) = cbind(pData(c17d), SVDF)
lkcd3 = gwSnpTests(probeId("GI_11038675-A")~X1+X2+X3+X4+
     X5+X6+X7+X8+X9+X10+X11+X12+X13+X14, c17d, chrnum("17"))
topSnps(lkcd3)


###################################################
### code chunk number 49: dotr (eval = FALSE)
###################################################
## library(GGtools)
## library(GGdata)
## t17f = transScores("GGdata", "17", rhs=~male, 
##   chrnames=c("chr1", "chr9"), wrapperEndo=
##   function(x) MAFfilter(x, low=.1), 
##   targdirpref="twfilt", geneApply=bestApply)
## intsave(t17f, file="t17f.rda")


###################################################
### code chunk number 50: gett (eval = FALSE)
###################################################
## topScores = function(tm) {
##  tm@base$scores[,1]/tm@base$shortfac
## }
## topGenes = function(tm) {
##  tm@base$guniv[ tm@base$inds[,1] ]
## }
## locusNames = function(tm) tm@base$snpnames
## geneNames = function(tm) tm@base$guniv
## geneIndcol = function(tm, col) tm@base$inds[,col]
## nthScores = function(tm, n) {
##  tm@base$scores[,n]/tm@base$shortfac
## }
## 


###################################################
### code chunk number 51: gettr
###################################################
tr17 = tr17_1_9()
tr17


###################################################
### code chunk number 52: getptr
###################################################
tr17_perm = tr17_1_9_perm()


###################################################
### code chunk number 53: lksco
###################################################
psco = topScores(tr17_perm)
summary(psco)
tp99 = quantile(psco, .99)


###################################################
### code chunk number 54: lktrrrr
###################################################
locw2 = which(nthScores(tr17,2)>tp99)
g1inds = geneIndcol(tr17,1)[locw2]
g2inds = geneIndcol(tr17,2)[locw2]
g1names = geneNames(tr17)[g1inds]
g2names = geneNames(tr17)[g2inds]
locnames = locusNames(tr17)[locw2]


###################################################
### code chunk number 55: lkm
###################################################
library(illuminaHumanv1.db)
nicevg = function(pr, rs, sms) {
 sym = get(pr, illuminaHumanv1SYMBOL)
 chr = get(pr, illuminaHumanv1CHR)
 plot_EvG( probeId(pr), rsid(rs), sms, 
   main=paste(sym, " chr", chr) )
}
par(mfrow=c(2,2))
nicevg( g1names[1], locnames[1], c17 )
nicevg( g2names[1], locnames[1], c17 )
nicevg( g1names[3], locnames[3], c17 )
nicevg( g2names[3], locnames[3], c17 )


###################################################
### code chunk number 56: getfdps
###################################################
pid = get("FDPS", revmap(illuminaHumanv1SYMBOL))
ind = which(geneNames(tr17) == pid)
isatop = any(geneIndcol(tr17,1) == ind)
isasec = any(geneIndcol(tr17,2) == ind)
isatop
isasec


###################################################
### code chunk number 57: makcplt (eval = FALSE)
###################################################
## c1 = getSS("cheung2010", "chr1")
## annotation(c1) = "hgfocus.db"
## annotation(CHE1) = "hgfocus.db"
## par(mfrow=c(2,2))
## plot_EvG(genesym("CHI3L2"), rsid("rs3934922"), CHE1)
## plot_EvG(genesym("CRYZ"), rsid("rs1475396"), CHE1)


###################################################
### code chunk number 58: dors
###################################################
library(Rsamtools)


###################################################
### code chunk number 59: doget
###################################################
bff=PileupFiles(dir(system.file("bam", package="ggtut"), 
   patt="bam$", full=TRUE))
# fix up names
no1 = gsub(".*bam.lit..", "", sapply(plpFiles(bff), path))
no2 = gsub("_1.bam", "", no1)
no3 = gsub(".bam", "", no2)
nanames = paste("NA", no3, sep="")
bff
# now have sample names


###################################################
### code chunk number 60: dobaca
###################################################
which = GRanges(seqnames="chr1", IRanges(111587452,width=1))
param = PileupParam(which=which)
pinfo = function(x) {
   tmp = apply(x[[3]], 2, function(z)
     z)
   tmp = tmp[c(2,3,5,9),]
   rownames(tmp)  = c("A", "C", "G", "T")
   tmp
   }
ans = applyPileups(bff, pinfo, param=param)[[1]]
colnames(ans) = nanames
ans


###################################################
### code chunk number 61: getchd
###################################################
c1 = getSS("cheung2010", "chr1")
c1s = smList(c1)[[1]]
rownames(c1s) = gsub("GM", "NA", c1$GMno)
c1ss = c1s[intersect(rownames(c1s),nanames),]
c1sss = c1ss[,"rs3934922"]
THESN = as(c1sss[,1], "character") 
rownames(THESN) = rownames(c1sss)
tt = as(THESN, "character")
names(tt) = rownames(c1ss)
ttt = tt[intersect(names(tt), colnames(ans))]
anss = ans[,names(ttt)]


###################################################
### code chunk number 62: lkh
###################################################
anssh = anss[,ttt=="A/B"]
plot(rep(1,ncol(anssh)),anssh[1,],xlim=c(.75,2.25),axes=FALSE,
xlab="allele rs3934922", ylab="CHI3L2 reads")
axis(1, at=c(1,2), labels=c("A", "C"))
axis(2)
segments(1,anssh[1,],2,anssh[2,])


###################################################
### code chunk number 63: lkd
###################################################
library(GenomicFeatures)
txdb = loadFeatures(system.file("sqlite/hg18.txdb.sqlite", 
 package="ggtut"))
allex = exons(txdb)


###################################################
### code chunk number 64: getse
###################################################
myr = GRanges(seqnames="chr1", IRanges(111500000, 111700000))


###################################################
### code chunk number 65: getex
###################################################
exinseg = allex[ allex %in% myr ]


###################################################
### code chunk number 66: dopar
###################################################
Pexinseg = PileupParam(which=exinseg)


###################################################
### code chunk number 67: findExWdiv
###################################################
hasTxDivers = function(x) {
   nuccts = x[[3]][c(2,3,5,9),,]
   apply(nuccts,3, function(w) any(apply(w,2,function(z)sum(z>5)>1)))
   }

ans = applyPileups(bff, hasTxDivers, param=Pexinseg)


###################################################
### code chunk number 68: eeee
###################################################
ediv = exinseg[ which(sapply(ans, any)) ]
ediv


###################################################
### code chunk number 69: getvarloc
###################################################
whichLocTxDivers = function(x) {
   nuccts = x[["seq"]][c(2,3,5,9),,]
   which(apply(apply(nuccts,3, function(w) apply(w,2,function(z)sum(z>5)>1)),2,any))
   }

Nparam = PileupParam(which=exinseg[which(sapply(ans,any))])

ans2 = applyPileups(bff, whichLocTxDivers, param=Nparam)


###################################################
### code chunk number 70: lklocs
###################################################
divlocs = lapply(1:length(ediv), function(x) (start(ediv[x]):end(ediv[x]))[ans2[[x]]])
divlocs
alldiv = unlist(divlocs)
length(alldiv)


###################################################
### code chunk number 71: lksnloc
###################################################
library(SNPlocs.Hsapiens.dbSNP.20090506)
c1locs = getSNPlocs("chr1")
sum(alldiv %in% c1locs$loc)
date()


###################################################
### code chunk number 72: getcalls (eval = FALSE)
###################################################
## library(GGtools)
## library(Rsamtools)
## exts  = seq(1, 80e6+10, by=10e6)
## st = exts[-length(exts)]
## en = exts[-1]-1
## # following file is 66 GB from 1000genomes.org
## tf = TabixFile("ALL.2of4intersection.20100804.genotypes.vcf.gz")
## gg = GRanges(seqnames="17", IRanges(st,en))
## for (i in 1:length(gg)) {
##   vv = vcf2sm(tf, gr=gg[i], nmetacol=9L)
##   intsave(vv, file=paste("vv", i, ".rda", sep=""))
## }


###################################################
### code chunk number 73: obtain (eval = FALSE)
###################################################
## getRules = function( tkgsm, basesm, locvec, try=200,
##   em.cntrl = c(1000,0.005, 1000, 0.005), use.hap = c(.99,.01) ) {
##  sntkg = colnames(tkgsm)
##  snbase = colnames(basesm)
## # unlocated loci need to be dropped
##  baseok = intersect(names(locvec),snbase)
##  tkok = intersect(names(locvec),sntkg)
##  tkgsm = tkgsm[, tkok]
##  basesm = basesm[, baseok]
##  sntkg = colnames(tkgsm)
##  snbase = colnames(basesm)
##  toimp = setdiff( sntkg, snbase )
##  usepred = setdiff( snbase, toimp )
##  yloc = locvec[ toimp ]
##  xloc = locvec[ usepred ]
##  rules = snp.imputation( basesm[,usepred], tkg[,toimp], xloc, yloc, try=200, 
##     em.cntrl = c(1000,0.005, 1000, 0.005), use.hap = c(.99,.01) )
## }
## 
## c17 = getSS("GGdata", "17", wrapperEndo=dropMonomorphies)
## base = smList(c17)[[1]]
## #
## # ceu1kg_17 was the restriction of the combined vv above
## # to 43 CEU individuals with 1000 genomes genotypes
## #
## tkg = get(load("ceu1kg_17.rda"))
## base = base[rownames(tkg),]
## 
## library(SNPlocs.Hsapiens.dbSNP.20090506)
## loc17 = getSNPlocs("chr17")
## snin = colnames(base)
## nams = paste("rs", loc17[,1], sep="")
## loc17 = loc17$loc
## names(loc17) = nams
## sninC = colnames(tkg)[ grep("chr", colnames(tkg)) ]
## kgloc = as.numeric(gsub("chr17:", "", sninC))
## names(kgloc) = sninC
## allloc = c(loc17, kgloc)
## 
## rules.n43 = getRules( tkg, base, allloc )