### R code from vignette source 'vignettes/GGtools/inst/doc/GGtools_2012.Rnw' ################################################### ### code chunk number 1: loadm ################################################### suppressPackageStartupMessages(library(GGtools)) library(parallel) ################################################### ### code chunk number 2: getd ################################################### g20 = getSS("GGdata", "20") ################################################### ### code chunk number 3: lkc ################################################### g20 class(g20) ################################################### ### code chunk number 4: doex ################################################### exprs(g20)[1:5,1:5] pData(g20)[1:4,] ################################################### ### code chunk number 5: lksm ################################################### smList(g20) as(smList(g20)[[1]][1:5,1:5], "matrix") ################################################### ### code chunk number 6: lksm2 ################################################### as(smList(g20)[[1]][1:5,1:5], "numeric") as(smList(g20)[[1]][1:5,1:5], "character") ################################################### ### code chunk number 7: dopeer (eval = FALSE) ################################################### ## library(peer) ## model = PEER() ## PEER_setPhenoMean(model, t(exprs(g20))) ## PEER_setNk(model, 10) ## PEER_setCovariates(model, matrix(1*g20$male,nc=1)) ## PEER_update(model) ## resid=t(PEER_getResiduals(model)) ## rownames(resid) = featureNames(g20) ## colnames(resid) = sampleNames(g20) ## g20peer10 = g20 ## g20peer10@assayData = assayDataNew("lockedEnvironment", exprs=resid) ################################################### ### code chunk number 8: dot ################################################### t1 = gwSnpTests(genesym("CPNE1")~male, g20, chrnum("20")) t1 topSnps(t1) ################################################### ### code chunk number 9: dopl ################################################### pdf(file="t1.pdf") plot(t1, snplocsDefault()) dev.off() pdf(file="t1evg.pdf") plot_EvG(genesym("CPNE1"), rsid("rs17093026"), g20) dev.off() ################################################### ### code chunk number 10: dol (eval = FALSE) ################################################### ## plot(t1, snplocsDefault()) ## plot_EvG(genesym("CPNE1"), rsid("rs17093026"), g20) ################################################### ### code chunk number 11: doloc (eval = FALSE) ################################################### ## library(snplocsDefault(), character.only=TRUE) ## S20 = getSNPlocs("ch20", as.GRanges=TRUE) ## GR20 = makeGRanges(t1, S20) ## library(rtracklayer) ## export(GR20, "~/cpne1new.wig") ################################################### ### code chunk number 12: prep ################################################### ################################################### ### code chunk number 13: do20 ################################################### g20 = GGtools:::restrictProbesToChrom(g20, "20") mads = apply(exprs(g20),1,mad) oo = order(mads, decreasing=TRUE) g20 = g20[oo[1:50],] tf = tempfile() dir.create(tf) e1 = eqtlTests(MAFfilter(g20, lower=0.05), ~male, geneApply=mclapply, targdir=tf) e1 ################################################### ### code chunk number 14: gettop ################################################### pm1 = colnames(e1@fffile) tops = sapply(pm1, function(x) topFeats(probeId(x), mgr=e1, n=1)) top6 = sort(tops, decreasing=TRUE)[1:6] ################################################### ### code chunk number 15: dopr6 ################################################### print(top6) ################################################### ### code chunk number 16: gettab ################################################### nms = strsplit(names(top6), "\\.") gn = sapply(nms,"[",1) sn = sapply(nms,"[",2) tab = data.frame(snp=sn,score=as.numeric(top6)) rownames(tab) = gn tab ################################################### ### code chunk number 17: ddstr ################################################### data(strMultPop) strMultPop[ strMultPop$rsid %in% tab$snp, ] ################################################### ### code chunk number 18: doc ################################################### if (file.exists("db2")) unlink("db2", recursive=TRUE) fn = probeId(featureNames(g20)) exTx = function(x) MAFfilter( x[fn, ], lower=0.05) b1 = best.cis.eQTLs("GGdata", ~male, radius=1e6, folderstem="db2", nperm=2, geneApply=mclapply, smFilter= exTx, chrnames="20", snpannopk=snplocsDefault()) ################################################### ### code chunk number 19: lkc ################################################### b1 ################################################### ### code chunk number 20: lkall ################################################### args(GGtools:::All.cis.eQTLs) ################################################### ### code chunk number 21: doall1 ################################################### b1.all = GGtools:::All.cis.eQTLs( maxfdr = 0.05, inbestcis = b1, smpack="GGdata", rhs=~1, chrnames = "20", smFilter4all = function(x) MAFfilter(x, lower=0.05)) ################################################### ### code chunk number 22: domo ################################################### if (file.exists("db2")) unlink("db2", recursive=TRUE) g20 = getSS("GGdata", "20") exTx = function(x) MAFfilter( clipPCs(x,1:10)[fn, ], lower=0.05) g20f = exTx(g20) ################################################### ### code chunk number 23: lkc2 ################################################### b2 = best.cis.eQTLs("GGdata", ~male, radius=1e6, folderstem="db2", nperm=2, geneApply=mclapply, smFilter= exTx, chrnames="20", snpannopk=snplocsDefault()) b2 ################################################### ### code chunk number 24: ggg ################################################### goodProbes = function(x) names(x@scoregr[elementMetadata(x@scoregr)$fdr<0.001]) ################################################### ### code chunk number 25: chkp ################################################### setdiff(goodProbes(b2), goodProbes(b1)) ################################################### ### code chunk number 26: domopic ################################################### newp = setdiff(goodProbes(b2), goodProbes(b1)) np = length(newp) bestSnp = function(pn, esm) elementMetadata(esm@scoregr[pn])$snpid par(mfrow=c(2,2)) plot_EvG(probeId(newp[1]), rsid(bestSnp(newp[1], b2)), g20, main="raw") plot_EvG(probeId(newp[1]), rsid(bestSnp(newp[1], b2)), g20f, main="PC-adjusted") plot_EvG(probeId(newp[np]), rsid(bestSnp(newp[np], b2)), g20, main="raw") plot_EvG(probeId(newp[np]), rsid(bestSnp(newp[np], b2)), g20f, main="PC-adjusted") ################################################### ### code chunk number 27: domoo ################################################### library(Biobase) suppressPackageStartupMessages(library(GGtools)) ex = matrix(0, nr=5, nc=3) pd = data.frame(v1 = 1:3, v2=5:7) colnames(ex) = rownames(pd) = LETTERS[1:3] adf = AnnotatedDataFrame(pd) rownames(ex) = letters[1:5] es = ExpressionSet(ex, phenoData=adf) exprs(es) pData(es) library(snpStats) mysnps = matrix(rep(1:3, 10), nr=3) # note 1=A/A, ... 0 = NA rownames(mysnps) = colnames(ex) mysm = new("SnpMatrix", mysnps) as(mysm, "character") as(mysm, "numeric") sml = make_smlSet(es, list(c1=mysm)) annotation(sml) colnames(smList(sml)[[1]]) ################################################### ### code chunk number 28: getss ################################################### toLatex(sessionInfo())