inst/doc/XDE.R

In XDE: XDE: a Bayesian hierarchical model for cross-study analysis of differential gene expression

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

###################################################
### code chunk number 1: XDE.Rnw:47-48
###################################################
options(width=60)


###################################################
### code chunk number 2: data
###################################################
library(XDE)
data(expressionSetList)
xlist <- expressionSetList
class(xlist)


###################################################
### code chunk number 3: validObject
###################################################
validObject(xlist)


###################################################
### code chunk number 4: covariate
###################################################
stopifnot(all(sapply(xlist, function(x, label){ label %in% varLabels(x)}, label="adenoVsquamous")))


###################################################
### code chunk number 5: initializeList
###################################################
params <- new("XdeParameter", esetList=xlist, phenotypeLabel="adenoVsquamous")
params


###################################################
### code chunk number 6: initializeList2 (eval = FALSE)
###################################################
## params <- new("XdeParameter", esetList=xlist,
## 	      phenotypeLabel="adenoVsquamous", one.delta=FALSE)


###################################################
### code chunk number 7: empiricalStart (eval = FALSE)
###################################################
## params <- new("XdeParameter", esetList=xlist, phenotypeLabel="adenoVsquamous", one.delta=FALSE)
## empirical <- empiricalStart(xlist, phenotypeLabel="adenoVsquamous")
## firstMcmc(params) <- empirical


###################################################
### code chunk number 8: burnin (eval = FALSE)
###################################################
## iterations(params) <- 3
## burnin(params) <- TRUE
## fit <- xde(params, xlist)


###################################################
### code chunk number 9: moreIterations (eval = FALSE)
###################################################
## iterations(params) <- 2
## fit2 <- xde(params, xlist, fit)


###################################################
### code chunk number 10: firstMcmc (eval = FALSE)
###################################################
## firstMcmc(params) <- lastMcmc(fit)
## seed(params) <- seed(fit)
## fit2 <- xde(params, xlist)


###################################################
### code chunk number 11: runMoreIterations (eval = FALSE)
###################################################
## burnin(params) <- FALSE
## iterations(params) <- 2000
## output(params)[c("potential", "acceptance",
## 		 "diffExpressed",
## 		 "nu", ##"DDelta",
##                  ##"delta",
## 		 "probDelta",
## 		 ##"sigma2",
## 		 "phi")] <- 0
## thin(params) <- 2
## directory(params) <- "logFiles"
## xmcmc <- xde(params, xlist)


###################################################
### code chunk number 12: savexmcmc (eval = FALSE)
###################################################
## ##this function should work
## postAvg <- calculatePosteriorAvg(xmcmc, NCONC=2, NDIFF=1)
## save(postAvg, file="~/madman/Rpacks/XDE/inst/logFiles/postAvg.rda")
## ##put browser in .standardizedDelta to fix tau2R, ...
## BES <- calculateBayesianEffectSize(xmcmc)
## save(BES, file="~/madman/Rpacks/XDE/inst/logFiles/BES.rda")
## save(xmcmc, file="~/madman/Rpacks/XDE/data/xmcmc.RData")
## q("no")


###################################################
### code chunk number 13: loadObject
###################################################
data(xmcmc, package="XDE")


###################################################
### code chunk number 14: savedParams
###################################################
output(xmcmc)[2:22][output(xmcmc)[2:22] == 1]


###################################################
### code chunk number 15: extractc2
###################################################
pathToLogFiles <- system.file("logFiles", package="XDE")
xmcmc@directory <- pathToLogFiles
c2 <- xmcmc$c2


###################################################
### code chunk number 16: c2
###################################################
par(las=1)
plot.ts(c2, ylab="c2", xlab="iterations", plot.type="single")


###################################################
### code chunk number 17: retrieveLogs
###################################################
getLogs <- function(object){
	params <- output(object)[output(object) == 1]
	params <- params[!(names(params) %in% c("nu", "phi", "DDelta", "delta", "sigma2", "diffExpressed"))]
	names(params)
}
param.names <- getLogs(xmcmc)
params <- lapply(lapply(as.list(param.names), function(name, object) eval(substitute(object$NAME_ARG, list(NAME_ARG=name))), object=xmcmc), as.ts)
names(params) <- param.names
tracefxn <- function(x, name)  plot(x, plot.type="single", col=1:ncol(x), ylab=name)
mapply(tracefxn, params, name=names(params))


###################################################
### code chunk number 18: bayesianEffectSize (eval = FALSE)
###################################################
## bayesianEffectSize(xmcmc) <- calculateBayesianEffectSize(xmcmc)


###################################################
### code chunk number 19: bes
###################################################
load(file.path(pathToLogFiles, "BES.rda"))
load(file.path(pathToLogFiles, "postAvg.rda"))


###################################################
### code chunk number 20: postAvgFig_conc
###################################################
par(las=1)
hist(postAvg[, "concordant"], breaks=50, xlim=c(0, 1), main="",
     xlab="posterior probability of concordant differential expression")


###################################################
### code chunk number 21: postAvgFig_disc
###################################################
par(las=1)
hist(postAvg[, "discordant"], breaks=50, xlim=c(0, 1), main="",
     xlab="posterior probability of discordant differential expression")


###################################################
### code chunk number 22: setPar
###################################################
op.conc <- symbolsInteresting(rankingStatistic=postAvg[, "concordant"])
op.disc <- symbolsInteresting(rankingStatistic=postAvg[, "discordant"])


###################################################
### code chunk number 23: conc
###################################################
par(las=1)
graphics:::pairs(BES[op.conc$order, ], pch=op.conc$pch, col=op.conc$col,
                 bg=op.conc$bg, upper.panel=NULL, cex=op.conc$cex)


###################################################
### code chunk number 24: disc
###################################################
graphics:::pairs(BES[op.disc$order, ], pch=op.disc$pch, col=op.disc$col, bg=op.disc$bg,
		 upper.panel=NULL, cex=op.disc$cex)


###################################################
### code chunk number 25: ssAlternatives (eval = FALSE)
###################################################
## ##t <- ssStatistic(statistic="t", phenotypeLabel="adenoVsquamous", esetList=xlist)
## tt <- rowttests(xlist, "adenoVsquamous", tstatOnly=TRUE)
## if(require(siggenes)){
##   sam <- ssStatistic(statistic="sam", phenotypeLabel="adenoVsquamous", esetList=xlist)
## }
## if(require(GeneMeta)){
##   z <- ssStatistic(statistic="z", phenotypeLabel="adenoVsquamous", esetList=xlist)
## }


###################################################
### code chunk number 26: tstatConc (eval = FALSE)
###################################################
## graphics:::pairs(tt[op.conc$order, ], pch=op.conc$pch, col=op.conc$col,
##                  bg=op.conc$bg, upper.panel=NULL, cex=op.conc$cex)
## 


###################################################
### code chunk number 27: zConc (eval = FALSE)
###################################################
## graphics:::pairs(tt[op.conc$order, ], pch=op.conc$pch, col=op.conc$col,
##                  bg=op.conc$bg, upper.panel=NULL, cex=op.conc$cex)
## 


###################################################
### code chunk number 28: tstatDisc (eval = FALSE)
###################################################
## graphics:::pairs(tt[op.disc$order, ], pch=op.disc$pch, col=op.disc$col,
##                  bg=op.disc$bg, upper.panel=NULL, cex=op.disc$cex)


###################################################
### code chunk number 29: pca (eval = FALSE)
###################################################
## tScores <- xsScores(tt, N=nSamples(xlist))
## samScores <- xsScores(sam, nSamples(xlist))
## zScores <- xsScores(z[, match(names(xlist), colnames(z))], N=nSamples(xlist))
## ##Concordant differential expression, we use the combined score from the random effects model directly
## zScores[, "concordant"] <- z[, "zSco"]


###################################################
### code chunk number 30: sessionInfo
###################################################
toLatex(sessionInfo())