###################################################
### chunk number 1: setup
###################################################
library("RbcBook1")
library("Biobase")
stopifnot(package.version("PROcess") >= package_version("1.3.4"))


###################################################
### chunk number 2: raw1
###################################################
library("PROcess")
fdat <- system.file("Test", package="PROcess")
fs <- list.files(fdat, pattern="\\.*csv\\.*", full.names=TRUE)
f1 <- read.files(fs[1])
plot(f1, type="l", xlab="m/z")
title(basename(fs[1]))


###################################################
### chunk number 3: bslnoffSingle
###################################################
bseoff <- bslnoff(f1,method="loess",bw=0.1,xlab="m/z",plot=TRUE)
title(basename(fs[1]))


###################################################
### chunk number 4: isPeak1
###################################################
pkgobj <- isPeak(bseoff,span=81,sm.span=11,
 plot=TRUE, zerothrsh=2,area.w = 0.003, ratio = 0.2, main="a)")


###################################################
### chunk number 5: isPeak2
###################################################
specZoom(pkgobj, xlim=c(5000,10000), main="b)")


###################################################
### chunk number 6: plotCalidef
###################################################
amu.cali <- c(1084,1638,3496,5807,7034)
plotCali <- function(f, main, lab.cali) {
  x <- read.files(f)
  plot(x, main=main, ylim=c(0, max(x[,2])), type="n")
  abline(h=0, col="gray")
  abline(v=amu.cali, col="salmon")
  if(lab.cali)
    axis(3, at=amu.cali,labels=amu.cali,las=3,
         tick=FALSE, col="salmon", cex.axis=0.94)
  lines(x)
  return(invisible(x))
}
dir.cali <- system.file("calibration", package="PROcess")
files    <- dir(dir.cali, full.names=TRUE)
i <- seq(along=files)


###################################################
### chunk number 7: plotCalido
###################################################
opar <- par(mfrow=c(4, 2), mar=c(2,2,3,1))
mapply(plotCali, files, LETTERS[i], i<=2) 
par(opar)


###################################################
### chunk number 8: plotCalishow eval=FALSE
###################################################
## mapply(plotCali, files, LETTERS[i], i<=2) 


###################################################
### chunk number 9: batchoff
###################################################
Mcal <- cache("Mcal", rmBaseline(dir.cali))


###################################################
### chunk number 10: batchoffshow eval=FALSE
###################################################
## Mcal <- rmBaseline(dir.cali)


###################################################
### chunk number 11: normalize
###################################################
M.r <- renorm(Mcal, cutoff = 400)


###################################################
### chunk number 12: cutoffDo
###################################################
cts <- round(10^(seq(2, 3.5, length=15)))
sdsFirst <- cache("sdsFirst", sapply(cts, avesd, Ma=Mcal))


###################################################
### chunk number 13: cutoffShow eval=FALSE
###################################################
## cts <- round(10^(seq(2, 4, length=14)))
## sdsFirst <- sapply(cts, avesd, Ma=Mcal)


###################################################
### chunk number 14: cutoff
###################################################
plot(cts, sdsFirst, xlab="cutpoint", pch=21, bg="red", log="x",
        ylab="average sd")


###################################################
### chunk number 15: getPeaksBatch.do
###################################################
peakfile <- "calipeak.csv"
if(!file.exists(peakfile))
   getPeaks(M.r, peakfile, ratio=0.1)


###################################################
### chunk number 16: getPeaksBatch.show eval=FALSE
###################################################
## peakfile <- "calipeak.csv"
## getPeaks(M.r, peakfile, ratio=0.1)


###################################################
### chunk number 17: QC
###################################################
qualRes <- quality(M.r, peakfile, cutoff=400)


###################################################
### chunk number 18: protobiomarkers
###################################################
bmkfile <- "calibmk.csv"
bmk1 <- pk2bmkr(peakfile,M.r,bmkfile, p.fltr=.5)
mk1 <- round(as.numeric(gsub("M","",names(bmk1))))
mk1


###################################################
### chunk number 19: plotCali2def
###################################################
plotCali2 <- function(...) {
  x <- plotCali(...)
  lines(x[,1]*2, x[,2]+25, col="blue")
}


###################################################
### chunk number 20: plotCali2do
###################################################
opar <- par(mfrow=c(4, 2), mar=c(2,2,3,1))
mapply(plotCali2, files, LETTERS[i], i<=2) 
par(opar)


###################################################
### chunk number 21: plotCali2show eval=FALSE
###################################################
## mapply(plotCali2, files, LETTERS[i], i<=2) 


###################################################
### chunk number 22: brcabseoff
###################################################
library("ProData")
f45c <- system.file("f45c", package="ProData")
fs <- dir(f45c,full.names=TRUE)


###################################################
### chunk number 23: rmBaseline.do
###################################################
M1 <- cache("M1", rmBaseline(f45c))


###################################################
### chunk number 24: rmBaseline.show eval=FALSE
###################################################
## M1 <- rmBaseline(f45c)


###################################################
### chunk number 25: brcacutoff
###################################################
data(f45cbmk)
SpecGrp <- pData(f45cbmk)
fns <- colnames(M1)
gi <-  regexpr("i+[0-9]+", fns)
specName <- substr(fns, gi, gi+attr(gi, "match.length")-1)
mt  <- match(SpecGrp[,2], toupper(specName))
M2 <- M1[, mt]
colnames(M2) <- SpecGrp[,2]


###################################################
### chunk number 26: sdsSecond.do
###################################################
stopifnot(!any(is.na(mt)))
sdsSecond <- cache("sdsSecond", 
  sapply(cts, avesd, Ma=M2[,SpecGrp[,1]=="D"]))


###################################################
### chunk number 27: sdsSecond.show eval=FALSE
###################################################
## sdsSecond <- sapply(cts, avesd, Ma=M2[,SpecGrp[,1]=="D"])


###################################################
### chunk number 28: brcacutoffFig
###################################################
plot(cts, sdsSecond, xlab="cutpoint", pch=21, bg="red", log="x",
        ylab="average sd")


###################################################
### chunk number 29: brcarenorm.cache
###################################################
nM <- cache("nM", renorm(M2, cutoff=1000))


###################################################
### chunk number 30: brcarenorm eval=FALSE
###################################################
## nM <- renorm(M2, cutoff=1000)


###################################################
### chunk number 31: brcapeakcallquality.do
###################################################
peakfile <- "f45cpeak.csv"
if(!file.exists(peakfile))
  getPeaks(nM, peakfile, ratio=.1)
qu <- cache("brcapeakcallquality", 
  quality(nM, peakfile, cutoff=1000))


###################################################
### chunk number 32: brcapeakcallquality.show eval=FALSE
###################################################
## peakfile <- "f45cpeak.csv"
## getPeaks(nM, peakfile, ratio=.1)
## qu <- quality(nM, peakfile, cutoff=1000)


###################################################
### chunk number 33: anybad1
###################################################
bad <- qu[,1] < 0.4 & qu[,2] <0.1 & qu[,3] < 1/2
sum(bad)


###################################################
### chunk number 34: anybad2
###################################################
stopifnot(!any(bad))


###################################################
### chunk number 35: brcanormalize
###################################################
drop <- SpecGrp[,1]=="D"
nM1 <- nM[,!drop]


###################################################
### chunk number 36: gelmap
###################################################
Ma <- binning(nM1, breaks=300)
colnames(Ma) <- SpecGrp[!drop,1]
par(xpd=TRUE)
marks <- c(2666, 5055, 7560, 7934)
sel   <- as.numeric(rownames(Ma))<10000
gelmap(log(Ma+1)[sel,],
        at.mz=marks, at.col=c(25, 90, 135),
        col=gray(10:0/10))
segments(x0=1013*c(1,1), y0=c(55,119),
         x1=10000*c(1,1), y1=c(55,119), col="red")
arrows(marks,c(-5,-5), marks, c(1,1),
        length=.08, angle=20, col="red")


###################################################
### chunk number 37: brcapeakalignDo
###################################################
peakfile1 <- "f45cpeak1.csv"
if( !file.exists(peakfile1) )
   getPeaks(nM1, peakfile1, ratio=.1)
bmkfile <- "f45cbmk.csv"
bmk <- cache("bmk", 
   pk2bmkr(peakfile1,nM1,bmkfile, p.fltr=0.1, eps=.003))


###################################################
### chunk number 38: brcapeakalignShow eval=FALSE
###################################################
## peakfile1 <- "f45cpeak1.csv"
## getPeaks(nM1, peakfile1, ratio=.1)
## bmkfile <- "f45cbmk.csv"
## bmk <- pk2bmkr(peakfile1,nM1,bmkfile, p.fltr=.1, eps=.003)


###################################################
### chunk number 39: brcabmk
###################################################
mks <- round(as.numeric(gsub("M","",names(bmk))))
length(mks)
mks[1:12]


###################################################
### chunk number 40: is.multiple
###################################################
mults <- is.multiple(mks, k=2:5)
## length(mults)
mults[[1]]