###################################################
### chunk number 1: 
###################################################
set.seed(1)
options(width=65)


###################################################
### chunk number 2: loadData
###################################################
library(les)
library(Biobase)
data(simTile)
treatment <- as.factor(phenoData(simTile)$condition == "treatment")
pos <- featureData(simTile)$position
exprs <- exprs(simTile)
regions <- c(100, 150, 600, 650)
cols <- rep(c("lightblue3", "lightgreen"), each=2)


###################################################
### chunk number 3: plotExpressionValues
###################################################
matplot(exprs, pch=".", xlab="Probe position", ylab="Expression")
abline(v=regions, col=cols)


###################################################
### chunk number 4: estimateProbeLevelStatistics
###################################################
library(limma)
design <- cbind(mean=1, diff=treatment)
fit <- lmFit(exprs, design)
fit <- eBayes(fit)
pval <- as.numeric(fit$p.value[, "diff"])


###################################################
### chunk number 5: plotProbeLevelStatistics
###################################################
plot(pos, pval, pch=".", xlab="Probe position", ylab="P-value")
abline(v=regions, col=cols)


###################################################
### chunk number 6: constructLes
###################################################
res <- Les(pos, pval)


###################################################
### chunk number 7: estimateLes
###################################################
win <- 30
res <- estimate(res, win)


###################################################
### chunk number 8: showPlotLes
###################################################
res
plot(res)
abline(v=regions, col=cols)


###################################################
### chunk number 9: estimateLes2
###################################################
win2 <- 50
res <- estimate(res, win2)


###################################################
### chunk number 10: showPlotLes2
###################################################
res
plot(res)
abline(v=regions, col=cols)


###################################################
### chunk number 11: threshold
###################################################
res <- threshold(res, grenander=TRUE, verbose=TRUE)


###################################################
### chunk number 12: regions
###################################################
theta <- 0.3
res <- regions(res, limit=theta, verbose=TRUE)
res
res["regions"]


###################################################
### chunk number 13: plotRegions
###################################################
region <- res["regions"]
borders <- c(region$start, region$end)
plot(res)
abline(v=regions, col=cols)
abline(v=borders, col="darkgray", lty=2)


###################################################
### chunk number 14: ci
###################################################
subset <- pos >= 580 & pos <= 670
res <- ci(res, subset, nBoot=50)


###################################################
### chunk number 15: plotCi
###################################################
plot(res, error="ci", limit=theta, xlim=c(580, 670))


###################################################
### chunk number 16: plotOptions
###################################################
plot(res, error="ci", region=TRUE, rug=TRUE, rugSide=3, main="LES for simulated data", probePch="*", probeCol="blue", errorCol="lightgray", regionCol="lightgreen", xlim=c(580, 670) , ylim=c(0, 0.8), limit=FALSE, lty=0)


###################################################
### chunk number 17: customWeightingFunction
###################################################
weightFoo <- function(distance, win) {
weight <- 1 - distance/win
return(weight)
}
res2 <- estimate(res, 20, weighting=weightFoo)


###################################################
### chunk number 18: sessionInfo
###################################################
toLatex(sessionInfo(), locale=FALSE)