## ----style, echo=FALSE, results='asis'-------------------------------------
BiocStyle::markdown()

## ----setup, echo=FALSE, message=FALSE--------------------------------------
library(Cardinal)
register(SerialParam())
options(Cardinal.verbose=FALSE)
options(Cardinal.progress=FALSE)
RNGkind("Mersenne-Twister")

## ----install, eval=FALSE---------------------------------------------------
#  install.packages("BiocManager")
#  BiocManager::install("Cardinal")

## ----library, eval=FALSE---------------------------------------------------
#  library(Cardinal)

## ----example-data----------------------------------------------------------
register(SerialParam())

set.seed(2020)
mse <- simulateImage(preset=1, npeaks=10, nruns=2, baseline=1)
mse

## ----pixel-data------------------------------------------------------------
pixelData(mse)

## ----coord-access----------------------------------------------------------
coord(mse)

## ----run-access------------------------------------------------------------
run(mse)[1:10]

## ----feature-data----------------------------------------------------------
featureData(mse)

## ----mz-access-------------------------------------------------------------
mz(mse)[1:10]

## ----image-data------------------------------------------------------------
imageData(mse)

## ----image-data-extract, eval=FALSE----------------------------------------
#  iData(mse, "intensity")

## ----spectra-access--------------------------------------------------------
spectra(mse)[1:5, 1:5]

## ----msi-continuous--------------------------------------------------------
mse
imageData(mse)

## ----msi-processed---------------------------------------------------------
set.seed(2020)
mse2 <- simulateImage(preset=3, npeaks=10, nruns=2)
mse3 <- as(mse2, "MSProcessedImagingExperiment")
mse3
imageData(mse3)

## ----resolution-access-----------------------------------------------------
resolution(mse3)

## ----resolution-update-1---------------------------------------------------
resolution(mse3) <- c(ppm=400)

## ----resolution-update-2---------------------------------------------------
dim(mse2)
dim(mse3)

## ----resolution-update-3---------------------------------------------------
mz(mse2)[1:10]
mz(mse3)[1:10]

## ----tiny-1----------------------------------------------------------------
set.seed(2020)
tiny <- simulateImage(preset=1, from=500, to=600, dim=c(3,3))
tiny

## ----tiny-2----------------------------------------------------------------
tiny2 <- as(tiny, "MSProcessedImagingExperiment")
tiny2

## ----writeMSIData-continous------------------------------------------------
path <- tempfile()
writeMSIData(tiny, file=path, outformat="imzML")

## ----writeMSIData-showfiles, echo=FALSE------------------------------------
grep(basename(path), list.files(dirname(path)), value=TRUE)

## ----writeMSIData-show-continuous-tag--------------------------------------
mzml <- readLines(paste0(path, ".imzML"))
grep("continuous", mzml, value=TRUE)

## ----writeMSIData-processed------------------------------------------------
path2 <- tempfile()
writeMSIData(tiny2, file=path2, outformat="imzML")

## ----writeMSIData-show-processed-tag---------------------------------------
mzml2 <- readLines(paste0(path2, ".imzML"))
grep("processed", mzml2, value=TRUE)

## ----tiny-3----------------------------------------------------------------
set.seed(2020)
tiny3 <- simulateImage(preset=1, from=500, to=600, dim=c(3,3), nruns=2)
runNames(tiny3)

## ----writeMSIData-multiple-runs--------------------------------------------
path3 <- tempfile()
writeMSIData(tiny3, file=path3, outformat="imzML")

## ----writeMSIData-multiple-runs-showfiles, echo=FALSE----------------------
grep(basename(path3), list.files(dirname(path3)), value=TRUE)

## ----readMSIData-continuos-------------------------------------------------
path_in <- paste0(path, ".imzML")
tiny_in <- readMSIData(path_in, attach.only=TRUE)
tiny_in

## ----readMSIData-spectra---------------------------------------------------
imageData(tiny_in)
spectra(tiny_in)

## ----readMSIData-spectra-2-------------------------------------------------
spectra(tiny_in)[1:5, 1:5]

## ----readMSIData-as-matrix-------------------------------------------------
spectra(tiny_in) <- as.matrix(spectra(tiny_in))
imageData(tiny_in)

## ----readMSIData-collect, eval=FALSE---------------------------------------
#  tiny_in <- collect(tiny_in)

## ----readMSIData-processed-------------------------------------------------
path2_in <- paste0(path2, ".imzML")
tiny2_in <- readMSIData(path2_in)
tiny2_in

## ----readMSIData-processed-2-----------------------------------------------
tiny2_in <- readMSIData(path2_in, mass.range=c(510,590),
						resolution=100, units="ppm")
tiny2_in

## ----readMSIData-processed-resolution--------------------------------------
resolution(tiny2_in) <- c(ppm=400)

## ----other-data------------------------------------------------------------
set.seed(2020)
s <- simulateSpectrum(n=9, peaks=10, from=500, to=600)

coord <- expand.grid(x=1:3, y=1:3)
run <- factor(rep("run0", nrow(coord)))

fdata <- MassDataFrame(mz=s$mz)
pdata <- PositionDataFrame(run=run, coord=coord)

out <- MSImagingExperiment(imageData=s$intensity,
							featureData=fdata,
							pixelData=pdata)
out

## ----plot-pixel, fig.height=3, fig.width=9---------------------------------
plot(mse, pixel=c(211, 611))

## ----plot-coord, fig.height=3, fig.width=9---------------------------------
plot(mse, coord=list(x=10, y=10), plusminus=1, fun=mean)

## ----plot-formula, fig.height=3, fig.width=9-------------------------------
plot(mse, intensity + I(-log1p(intensity)) ~ mz, pixel=211, superpose=TRUE)

## ----image-mz, fig.height=4, fig.width=9-----------------------------------
image(mse, mz=1200)

## ----image-plusminus, fig.height=4, fig.width=9----------------------------
image(mse, mz=1227, plusminus=0.5, fun=mean)

## ----image-run, fig.height=4, fig.width=5----------------------------------
image(mse, mz=1227, subset=run(mse) == "run0")

## ----image-subset, fig.height=8, fig.width=9-------------------------------
image(mse, mz=c(1200, 1227), subset=circle)

## ----image-smooth, fig.height=4, fig.width=9-------------------------------
image(mse, mz=1200, smooth.image="gaussian")

## ----image-contrast, fig.height=4, fig.width=9-----------------------------
image(mse, mz=1200, contrast.enhance="histogram")

## ----image-colorscale, fig.height=4, fig.width=9---------------------------
image(mse2, mz=1136, colorscale=magma)

## ----image-layout, fig.height=2, fig.width=9-------------------------------
image(mse2, mz=c(781, 1361), layout=c(1,4), colorscale=magma)

## ----image-superpose, fig.height=4, fig.width=9----------------------------
image(mse2, mz=c(781, 1361), superpose=TRUE, normalize.image="linear")

## ----image-formula, fig.height=4, fig.width=9------------------------------
image(mse2, log1p(intensity) ~ x * y, mz=1136, colorscale=magma)

## ----image3d---------------------------------------------------------------
set.seed(1)
mse3d <- simulateImage(preset=9, nruns=7, dim=c(7,7), npeaks=10,
						peakheight=c(2,4), representation="centroid")

image3D(mse3d, mz=c(1001, 1159), colorscale=plasma, cex=3, theta=30, phi=30)

## ----select-ROI, eval=FALSE------------------------------------------------
#  sampleA <- selectROI(mse, mz=1200, subset=run(mse) == "run0")
#  sampleB <- selectROI(mse, mz=1200, subset=run(mse) == "run1")

## ----makeFactor, eval=FALSE------------------------------------------------
#  regions <- makeFactor(A=sampleA, B=sampleB)

## ----pdf, eval=FALSE-------------------------------------------------------
#  pdf_file <- paste0(tempfile(), ".pdf")
#  
#  pdf(file=pdf_file, width=9, height=4)
#  image(mse, mz=1200)
#  dev.off()

## ----dark-mode-1, fig.height=4, fig.width=9--------------------------------
darkmode()
image(mse, mz=1200)

## ----dark-mode-2, fig.height=4, fig.width=9--------------------------------
darkmode()
image(mse2, mz=1135, colorscale=magma)

## ----light-mode------------------------------------------------------------
lightmode()

## ----print, eval=FALSE-----------------------------------------------------
#  g <- image(mse, mz=1200)
#  print(g)

## ----subset-1--------------------------------------------------------------
# subset first 10 mass features and first 5 pixels
mse[1:10, 1:5]

## ----features--------------------------------------------------------------
# get row index corresponding to m/z 1200
features(mse, mz=1200)

# get row indices corresponding to m/z 1199 - 1201
features(mse, 1199 < mz & mz < 1201)

## ----pixels----------------------------------------------------------------
# get column indices corresponding to x = 10, y = 10 in all runs
pixels(mse, coord=list(x=10, y=10))

# get column indices corresponding to x <= 3, y <= 3 in "run0"
pixels(mse, x <= 3, y <= 3, run == "run0")

## ----subset-2--------------------------------------------------------------
fid <- features(mse, 1199 < mz & mz < 1201)
pid <- pixels(mse, x <= 3, y <= 3, run == "run0")
mse[fid, pid]

## ----slice-----------------------------------------------------------------
# slice image for first mass feature
a <- slice(mse, 1)
dim(a)

## ----slice-mz--------------------------------------------------------------
# slice image for m/z 1200
a2 <- slice(mse, mz=1200, .preserve=TRUE)
dim(a2)

## ----slice-image, fig.height=4, fig.width=4--------------------------------
image(a2[,,1,1], col=bw.colors(100))

## ----cbind-divide----------------------------------------------------------
# divide dataset into separate objects for each run
mse_run0 <- mse[,run(mse) == "run0"]
mse_run1 <- mse[,run(mse) == "run1"]
mse_run0
mse_run1

## ----cbind-recombine-------------------------------------------------------
# recombine the separate datasets back together
mse_cbind <- cbind(mse_run0, mse_run1)
mse_cbind

## ----pData-set-------------------------------------------------------------
mse$region <- makeFactor(circle=mse$circle,
							bg=!mse$circle)
pData(mse)

## ----iData-set-------------------------------------------------------------
iData(mse, "log2intensity") <- log2(iData(mse) + 1)
imageData(mse)

## ----spectra-get-----------------------------------------------------------
spectra(mse, "log2intensity")[1:5, 1:5]

## ----centroided-get--------------------------------------------------------
centroided(mse)

## ----centroided-set, eval=FALSE--------------------------------------------
#  centroided(mse) <- FALSE

## ----dplyr, eval=FALSE-----------------------------------------------------
#  # filter mass range
#  filter(mse, mz > 700, mz < 900)
#  
#  # select pixels
#  select(mse, x <= 3, y <= 3, run == "run0")
#  
#  # chain verbs together
#  mse %>%
#  	filter(mz < 1000) %>%
#  	select(x == 10, y == 10)
#  
#  # calculate mean spectrum
#  summarize(mse, .stat="mean", .by="feature", .as="DataFrame")
#  
#  # calculate tic
#  summarize(mse, .stat=c(tic="sum"), .by="pixel", .as="DataFrame")
#  
#  # calculate mean spectrum of circle region
#  mse %>%
#  	select(region == "circle") %>%
#  	summarize(.stat="mean", .as="DataFrame")

## ----matter----------------------------------------------------------------
# coerce spectra to file-based matter matrix
spectra(mse) <- matter::as.matter(spectra(mse))

spectra(mse)
imageData(mse)

## --------------------------------------------------------------------------
mse <- collect(mse)
imageData(mse)

## ----eval=FALSE------------------------------------------------------------
#  mse3 <- collect(mse3, as.matrix=TRUE)

## ----coerce----------------------------------------------------------------
mse3
as(mse3, "MSContinuousImagingExperiment")

## ----coerce-2, eval=FALSE--------------------------------------------------
#  # requires CardinalWorkflows installed
#  data(cardinal, package="CardinalWorkflows")
#  cardinal2 <- as(cardinal, "MSImagingExperiment")

## ----process-1-------------------------------------------------------------
mse_tic <- process(mse, function(x) length(x) * x / sum(x), label="norm")
mse_tic

## ----process-2-------------------------------------------------------------
mse_pre <- mse %>%
	process(function(x) length(x) * x / sum(x), label="norm", delay=TRUE) %>%
	process(function(x) signal::sgolayfilt(x), label="smooth", delay=TRUE)

processingData(mse_pre)

## ----process-3-------------------------------------------------------------
mcols(processingData(mse_pre))[,-1]

## ----process-4-------------------------------------------------------------
mse_proc <- process(mse_pre)
mse_proc

## ----process-5-------------------------------------------------------------
mse_pre %>%
	select(x <= 3, y <= 3) %>%
	filter(mz <= 1000) %>%
	process()

## ----normalize-------------------------------------------------------------
mse_pre <- normalize(mse, method="tic")

## ----smoothSignal-plot, fig.height=7, fig.width=9, results='hide'----------
mse %>% smoothSignal(method="gaussian") %>%
	select(x==10, y==10, run=="run0") %>%
	process(plot=TRUE,
			par=list(main="Gaussian smoothing", layout=c(3,1)))

mse %>% smoothSignal(method="ma") %>%
	select(x==10, y==10, run=="run0") %>%
	process(plot=TRUE, par=list(main="Moving average smoothing"))

mse %>% smoothSignal(method="sgolay") %>%
	select(x==10, y==10, run=="run0") %>%
	process(plot=TRUE, par=list(main="Savitzky-Golay smoothing"))

## ----smoothSignal----------------------------------------------------------
mse_pre <- smoothSignal(mse_pre, method="gaussian")

## ----reduceBaseline-plot, fig.height=5, fig.width=9, results='hide'--------
mse %>% reduceBaseline(method="locmin") %>%
	select(x==10, y==10, run=="run0") %>%
	process(plot=TRUE, par=list(main="Local minima", layout=c(2,1)))

mse %>% reduceBaseline(method="median") %>%
	select(x==10, y==10, run=="run0") %>%
	process(plot=TRUE, par=list(main="Median interpolation"))

## ----reduceBaseline--------------------------------------------------------
mse_pre <- reduceBaseline(mse_pre, method="locmin")

## ----process-mse, eval=FALSE-----------------------------------------------
#  mse_pre <- process(mse_pre)

## ----peakPick-plot, fig.height=7, fig.width=9, results='hide'--------------
mse_pre %>% select(x==10, y==10, run=="run0") %>%
	process() %>%
	peakPick(method="mad") %>%
	process(plot=TRUE, par=list(main="MAD noise", layout=c(3,1)))

mse_pre %>% select(x==10, y==10, run=="run0") %>%
	process() %>%
	peakPick(method="simple") %>%
	process(plot=TRUE,
			par=list(main="Simple SD noise"))

mse_pre %>% select(x==10, y==10, run=="run0") %>%
	process() %>%
	peakPick(method="adaptive") %>%
	process(plot=TRUE, par=list(main="Adaptive SD noise"))

## ----peakPick--------------------------------------------------------------
mse_peaks <- peakPick(mse_pre, method="mad")

## ----peakAlign-------------------------------------------------------------
mse_peaks <- peakAlign(mse_pre, tolerance=200, units="ppm")

## ----peakFilter------------------------------------------------------------
mse_peaks <- peakFilter(mse_pre, freq.min=0.05)

## ----peakBin, eval=FALSE---------------------------------------------------
#  mse_peaks <- process(mse_peaks)
#  mse_peaks <- peakBin(mse_pre, ref=mz(mse_peaks), type="height")
#  mse_peaks <- mse_peaks %>% process()

## ----peakBin-2, fig.height=3, fig.width=9----------------------------------
mzref <- mz(metadata(mse)$design$featureData)
mse_peaks <- peakBin(mse_pre, ref=mzref, type="height")
mse_peaks <- mse_peaks %>% select(x %in% 9:11, y %in% 9:11) %>% process()
plot(mse_peaks, coord=list(x=10, y=10))

## ----unaligned-spectra, fig.height=3, fig.width=9--------------------------
set.seed(2020)
mse4 <- simulateImage(preset=1, npeaks=10, from=500, to=600,
							sdmz=750, units="ppm")

plot(mse4, pixel=185:195, xlim=c(535, 570), key=FALSE, superpose=TRUE)

## ----mzAlign, fig.height=3, fig.width=9------------------------------------
mse4_align <- mzAlign(mse4, tolerance=2000, units="ppm")
mse4_align <- process(mse4_align)
plot(mse4_align, pixel=185:195, xlim=c(535, 570), key=FALSE, superpose=TRUE)

## ----mzBin, fig.height=3, fig.width=9--------------------------------------
mse_bin <- mzBin(mse_pre, from=1000, to=1600, resolution=1000, units="ppm")
mse_bin <- select(mse_bin, x %in% 9:11, y %in% 9:11) %>% process()
plot(mse_bin, coord=list(x=10, y=10))

## ----mzFilter, fig.height=3, fig.width=9-----------------------------------
mse_filt <- mzFilter(mse_pre, thresh.max=0.05)
mse_filt <- select(mse_filt, x %in% 9:11, y %in% 9:11) %>% process()
plot(mse_filt, coord=list(x=10, y=10), type='h')

## ----process-workflow, fig.height=5, fig.width=9, results='hide'-----------
mse_proc <- mse %>%
	normalize() %>%
	smoothSignal() %>%
	reduceBaseline() %>%
	peakPick()

# preview processing
mse_proc %>%
	select(x == 10, y == 10, run == "run0") %>%
	process(plot=TRUE, par=list(layout=c(2,2)))

## ----process-workflow-2, eval=FALSE----------------------------------------
#  # process detected peaks
#  mse_peakref <- mse_proc %>%
#  	peakAlign() %>%
#  	peakFilter() %>%
#  	process()
#  
#  # bin profile spectra to peaks
#  mse_peaks <- mse %>%
#  	normalize() %>%
#  	smoothSignal() %>%
#  	reduceBaseline() %>%
#  	peakBin(mz(mse_peakref))

## ----pixelApply, fig.height=4, fig.width=9---------------------------------
# calculate the TIC for every pixel
tic <- pixelApply(mse, sum)

image(mse, tic ~ x * y)

## ----featureApply, fig.height=3, fig.width=9-------------------------------
# calculate the mean spectrum
smean <- featureApply(mse, mean)

plot(mse, smean ~ mz)

## ----spatialApply, fig.height=4, fig.width=9-------------------------------
# calculate a spatially-smoothed TIC
sptic <- spatialApply(mse, .r=1, .fun=mean)

image(mse, sptic ~ x * y)

## ----eval=FALSE------------------------------------------------------------
#  # run serially, rather than in parallel
#  tic <- pixelApply(mse, sum, BPPARAM=SerialParam())

## ----registered------------------------------------------------------------
registered()

## ----register, eval=FALSE--------------------------------------------------
#  # register a SNOW backend
#  register(SnowParam())

## ----RNGkind, eval=FALSE---------------------------------------------------
#  RNGkind("L'Ecuyer-CMRG")
#  set.seed(1)

## ----session-info----------------------------------------------------------
sessionInfo()