### R code from vignette source 'gsri.Rnw' ### Encoding: UTF-8 ################################################### ### code chunk number 1: settings ################################################### set.seed(1) options(width=65, SweaveHooks=list(fig=function() par(mar=c(5.1, 5.1, 4.1, 1.1)))) ################################################### ### code chunk number 2: preparePdfCdfFig ################################################### x <- seq(0, 1, len=50) r <- 0.2 rate <- 10 d0 <- dunif(x) d1 <- dexp(x, rate) d <- (1-r)*d0 + r*d1 c0 <- punif(x) c1 <- pexp(x, rate) c <- (1-r)*c0 + r*c1 x <- seq(0, 1, len=50) cex <- 1.5 ################################################### ### code chunk number 3: pdfFig ################################################### getOption("SweaveHooks")[["fig"]]() plot(x, d, type="n", xaxt="n", yaxt="n", ylim=c(0, max(d)), xlab=expression(paste("p-value ", italic(p))), ylab=expression(paste("probability density ", rho(p))), cex.lab=cex) lines(c(0, 1), rep(1-r, 2), lwd=2, col="darkgray") lines(x, d, lwd=2) axis(1, at=seq(0, 1, len=5), labels=c(0, NA, 0.5, NA, 1), cex.axis=cex) axis(2, at=seq(0, max(d), by=0.5), labels=c(0, NA, 1, NA, 2, NA), cex.axis=cex) axis(2, at=1-r, labels=expression(paste(1-italic(r))), cex.axis=cex) text(0.5, (1-r)/2, labels=expression(1-italic(r)), cex=cex, offset=NULL, adj=c(0.5, 0.5)) text(0.09, 1.05, labels=expression(italic(r)), cex=cex, offset=NULL, adj=c(0, 0)) ################################################### ### code chunk number 4: cdfFig ################################################### getOption("SweaveHooks")[["fig"]]() plot(x, c, type="n", xaxt="n", yaxt="n", xlim=c(0, 1), ylim=c(0, 1), xlab=expression(paste("p-value ", italic(p))), ylab=expression(paste("cumulative density ", F(p))), cex.lab=cex) lines(c(0, 1), c(r, 1), lwd=2, col="darkgray") lines(x, c, lwd=2) axis(1, at=seq(0, 1, len=5), labels=c(0, NA, 0.5, NA, 1), cex.axis=cex) axis(2, at=seq(0, 1, len=5), labels=c(0, NA, 0.5, NA, 1), cex.axis=cex) axis(2, at=r, labels=expression(paste(italic(r))), cex.axis=cex) text(0.75, 0.75, labels=expression(1-italic(r)), cex=cex, offset=NULL, adj=c(0, 0)) ################################################### ### code chunk number 5: extractData ################################################### library(Biobase) data(sample.ExpressionSet) eset <- sample.ExpressionSet eset exprs <- exprs(eset) phenotypes <- pData(phenoData(eset)) summary(phenotypes) ################################################### ### code chunk number 6: gsriAllProbes ################################################### library(GSRI) gAllProbes <- gsri(eset, phenotypes$type) gAllProbes ################################################### ### code chunk number 7: gsriAllGenesSet ################################################### library(GSEABase) gs <- GeneSet(eset, setName="allGenes") ind <- grep("^AFFX", geneIds(gs), invert=TRUE) gsAllGenes <- gs[ind] gsAllGenes ################################################### ### code chunk number 8: gsriAllGenes ################################################### gAllGenesType <- gsri(eset, phenotypes$type, gsAllGenes, name="allGenesType") gAllGenesSex <- gsri(exprs, phenotypes$sex, gsAllGenes, name="allGenesSex") gAllGenesType gAllGenesSex ################################################### ### code chunk number 9: gsriChr17 ################################################### gsChr17 <- getBroadSets(system.file("extdata", "c1chr17.xml", package="GSRI")) gsChr17 gChr17 <- gsri(eset, phenotypes$type, gsChr17) gChr17 ################################################### ### code chunk number 10: gsriCol5 ################################################### gmt <- getGmt(system.file("extdata", "c1c10.gmt", package="GSRI")) gCol5 <- gsri(eset, phenotypes$type, gmt) gCol5 gCol5Sort <- sortGsri(gCol5, c("nRegGenes", "pRegGenes")) summary(gCol5Sort) exportFile <- tempfile() export(gCol5Sort, exportFile) ################################################### ### code chunk number 11: gsriScoreFtest ################################################### phenotypes$class <- cut(phenotypes$score, seq(0, 1, length.out=4), label=c("low", "medium", "high")) summary(phenotypes$class) g3 <- gsri(eset, phenotypes$class, gsChr17, test=rowF) g3 ################################################### ### code chunk number 12: gsriAllGenesArgs ################################################### g3arg2 <- gsri(eset, phenotypes$class, gsChr17, test=rowF, name="chr17_2", nBoot=50, alpha=0.1, grenander=FALSE) g3arg2 ################################################### ### code chunk number 13: limmaTest ################################################### library(limma) limmaTest <- function(exprs, groups, id, index, testArgs) { design <- cbind(offset=1, diff=groups) fit <- lmFit(exprs[ ,index], design) fit <- eBayes(fit) pval <- fit$p.value[id,"diff"] return(pval) } g3Limma <- gsri(eset, phenotypes$type, gsChr17, test=limmaTest) g3Limma ################################################### ### code chunk number 14: plot1 ################################################### getOption("SweaveHooks")[["fig"]]() plot(g3) ################################################### ### code chunk number 15: plot2 ################################################### getOption("SweaveHooks")[["fig"]]() plot(gCol5, 5, ecdf=list(type="o"), plot=list(xlab="p", main="GSRI results: chr9"), reg=list(col="lightblue", lty=1, lwd=1.5), gsri=list(col="darkblue")) ################################################### ### code chunk number 16: weights ################################################### library(hgu95av2.db) gNames <- rownames(exprs(eset)) ind <- Lkeys(hgu95av2GO) %in% gNames evidence <- factor(toTable(hgu95av2GO)[ind,"Evidence"]) summary(evidence) l <- lapply(gNames, function(name, names, evidence) evidence[names %in% name], gNames, evidence) expInd <- sapply(l, function(l) any(l %in% "EXP")) goWeight <- rep(0.5, length.out=length(expInd)) goWeight[expInd] <- 1 gCol5go <- gsri(eset, phenotypes$type, weight=goWeight) gCol5go gCol5go2 <- gsri(eset, phenotypes$type, gmt, weight=goWeight) gCol5go2 ################################################### ### code chunk number 17: sessionInfo ################################################### toLatex(sessionInfo(), locale=FALSE)