### R code from vignette source 'poster.Rnw'
### Encoding: UTF-8

###################################################
### code chunk number 1: style-Sweave
###################################################
BiocStyle::latex_old()
options(prompt = "R> ", continue = "+  ", width = 60, useFancyQuotes = FALSE)


###################################################
### code chunk number 2: poster.Rnw:130-150 (eval = FALSE)
###################################################
## library(cosmiq)
## cdfpath <- file.path(find.package("faahKO"),
## 	"cdf")
## 
## my.input.files <- dir(c(paste(cdfpath,
## 	"WT", sep='/'), 
## 	paste(cdfpath, "KO", sep='/')),
## 	full.names=TRUE)
## 
## # run cosmiq wrapper function
## #
## x <- cosmiq(files = my.input.files, 
## 	mzbin=0.25, 
## 	SNR.Th=0, 
## 	linear=TRUE)
## 
## # graph result
## image(t(x$eicmatrix), 
## 	main='mz versus RT map')
## head(x$xs@peaks)


###################################################
### code chunk number 3: poster.Rnw:171-172
###################################################
options(width=80)


###################################################
### code chunk number 4: poster.Rnw:181-194
###################################################
library(cosmiq)
cdfpath <- file.path(find.package("faahKO"),
	"cdf")
my.input.files <- dir(c(paste(cdfpath,
	"WT", sep='/'), 
	paste(cdfpath, "KO", sep='/')),
	full.names=TRUE)

#
# create xcmsSet object
# todo
xs <- new("xcmsSet")
xs@filepaths <- my.input.files


###################################################
### code chunk number 5: poster.Rnw:200-204
###################################################
class <- as.data.frame(c(rep("KO",6),
	rep("WT", 6)))
rownames(class) <- basename(my.input.files)
xs@phenoData <- class


###################################################
### code chunk number 6: poster.Rnw:210-212
###################################################
xs.attr <- attributes(xs)
xs.attr$phenoData


###################################################
### code chunk number 7: combined
###################################################

x <- combine_spectra(xs=xs, mzbin=0.25, 
	linear=TRUE, continuum=FALSE)

plot(x$mz, x$intensity, type='l',
	main='combined spectra',
	xlab='m/Z', ylab='ion intensity')



###################################################
### code chunk number 8: plot2
###################################################

xy <- peakdetection(x=x$mz, y=x$intensity, 
    scales=1:10, 
    SNR.Th=1.0, 
    SNR.area=20, mintr=0.5)

id.peakcenter<-xy[,4]

filter.mz <- 400 < x$mz & x$mz < 450
plot(x$mz[filter.mz], 
    x$intensity[filter.mz], 
    main='Detection of relevant masses',
    type='l',
    xlab='m/Z', 
    ylab='ion intensity')

points(x$mz[id.peakcenter], 
    x$intensity[id.peakcenter], 
    col='red', type='h')



###################################################
### code chunk number 9: poster.Rnw:301-315
###################################################
# create dummy object
xs@peaks <- matrix(c(rep(1,
	length(my.input.files) * 6), 
	1:length(my.input.files)),
	ncol=7)

colnames(xs@peaks) <- c("mz",
	"mzmin", "mzmax", "rt", 
	"rtmin", "rtmax", "sample")

xs <- xcms::retcor(xs,
	method = "obiwarp", profStep=1, 
	distFunc="cor", center=1)



###################################################
### code chunk number 10: image
###################################################
eicmat <- 
eicmatrix(xs=xs, xy=xy, center=1)

# process a reduced mz range 
# for a better package build performance
(eicmat.mz.range 
<- range(which(475 < xy[,1] 
	& xy[,1] < 485)))

eicmat.filter <-
eicmat[eicmat.mz.range[1]:
eicmat.mz.range[2], ]

xy.filter <- xy[eicmat.mz.range[1]:eicmat.mz.range[2], ]

# determine the new range and 
# plot the mz versus RT map
(rt.range <- 
range(as.double(colnames(eicmat.filter))))
(mz.range <-
range(as.double(row.names(eicmat.filter))))

image(log(t(eicmat.filter))/log(2), 
    main = 'overlay of 12
    samples using faahKO',
    col = rev(gray(1:20/20)),
    xlab = 'rt [in seconds]', 
    ylab = 'm/z', 
    axes = FALSE)

axis(1, seq(0,1, length=6), 
    round(seq(rt.range[1], 
    rt.range[2], 
    length = 6)))

axis(2, seq(0,1, length = 4), 
    round(seq(mz.range[1], 
    mz.range[2], 
    length = 4), 2))

#
# determine the chromatographic peaks
rxy <- retention_time(xs=xs, 
    RTscales = c(1:10, seq(12,32, by = 2)), 
    xy = xy.filter, 
    eicmatrix = eicmat.filter, 
    RTSNR.Th = 120, 
    RTSNR.area = 20)

rxy.rt <- 
(rxy[,4] - rt.range[1]) / diff(rt.range)

rxy.mz <- 
(rxy[,1] - mz.range[1]) / diff(mz.range)

points(rxy.rt, rxy.mz, 
	pch = 16, 
	lwd = 2, 
	col = rgb(0.9, 0.1, 0.1, alpha=0.5))


###################################################
### code chunk number 11: poster.Rnw:410-412
###################################################
xs <- create_datamatrix(xs = xs,
	rxy = rxy)


###################################################
### code chunk number 12: poster.Rnw:422-423
###################################################
peaktable <- xcms::peakTable(xs)


###################################################
### code chunk number 13: poster.Rnw:426-427
###################################################
idx <- order(rowSums(peaktable[,8:19]), decreasing=TRUE)


###################################################
### code chunk number 14: xtable1
###################################################
library(xtable)
xtab <- peaktable[idx[1:20], ]
print.xtable(xtable(xtab, 
  caption="The spreadsheet shows the top 20 most intense rows (\\Rcode{order(rowSums(peaktable[,8:19]), decreasing=TRUE)}) of the \\Rcode{peakTable} result.", 
	label="Table:xtable1", 
	include.rownames=FALSE), 
	table.placement = "H",
	scalebox="0.50")