inst/doc/Cardinal-2-guide.R

In Cardinal: A mass spectrometry imaging toolbox for statistical analysis

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

## ----setup, echo=FALSE, message=FALSE-----------------------------------------
library(Cardinal)
setCardinalBPPARAM(SerialParam())
setCardinalVerbose(FALSE)
RNGkind("Mersenne-Twister")

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

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

## ----example-data-------------------------------------------------------------
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, drop=FALSE)
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

## ----manip--------------------------------------------------------------------
# subset by mass range
subsetFeatures(mse, mz > 700, mz < 900)

# subset by pixel coordinates
subsetPixels(mse, x <= 3, y <= 3, run == "run0")

# subset by mass range + pixel coordinates
subset(mse, mz > 700 & mz < 900, x <=3 & y <= 3 & run == "run0")

# chain verbs together
mse %>%
	subsetFeatures(mz > 700, mz < 900) %>%
	subsetPixels(x <= 3, y <= 3, run == "run0")

# calculate mean spectrum
summarizeFeatures(mse, "mean", as="DataFrame")

# calculate tic
summarizePixels(mse, c(tic="sum"), as="DataFrame")

# calculate mean spectrum of circle region
mse %>%
	subsetPixels(region == "circle") %>%
	summarizeFeatures("mean", as="DataFrame")

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

spectra(mse)
imageData(mse)

## -----------------------------------------------------------------------------
mse <- pull(mse)
imageData(mse)

## ----eval=FALSE---------------------------------------------------------------
#  mse3 <- pull(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 %>%
	subsetPixels(x <= 3, y <= 3) %>%
	subsetFeatures(mz <= 1000) %>%
	process()

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

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

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

mse %>% smoothSignal(method="sgolay") %>%
	subsetPixels(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") %>%
	subsetPixels(x==10, y==10, run=="run0") %>%
	process(plot=TRUE, par=list(main="Local minima", layout=c(2,1)))

mse %>% reduceBaseline(method="median") %>%
	subsetPixels(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 %>% subsetPixels(x==10, y==10, run=="run0") %>%
	process() %>%
	peakPick(method="mad") %>%
	process(plot=TRUE, par=list(main="MAD noise", layout=c(3,1)))

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

mse_pre %>% subsetPixels(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 %>% subsetPixels(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)

## ----mzAlign1, fig.height=3, fig.width=9--------------------------------------
mse4_mean <- summarizeFeatures(mse4)
mse4_peaks <- peakPick(mse4_mean, SNR=2)
mse4_peaks <- peakAlign(mse4_peaks, tolerance=1000, units="ppm")
mse4_peaks <- process(mse4_peaks)
fData(mse4_peaks)

mse4_align1 <- mzAlign(mse4, ref=mz(mse4_peaks), tolerance=2000, units="ppm")
mse4_align1 <- process(mse4_align1)
plot(mse4_align1, pixel=185:195, xlim=c(535, 570), key=FALSE, superpose=TRUE)

## ----mzAlign2, fig.height=3, fig.width=9--------------------------------------
mse4_align2 <- mzAlign(mse4, tolerance=2000, units="ppm")
mse4_align2 <- process(mse4_align2)
plot(mse4_align2, 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 <- subsetPixels(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, freq.min=0.05)
mse_filt <- subsetPixels(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 %>%
	subsetPixels(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 in parallel, rather than serially
#  tic <- pixelApply(mse, sum, BPPARAM=MulticoreParam())

## ----registered---------------------------------------------------------------
BiocParallel::registered()

## ----getCardinalBPPARAM-------------------------------------------------------
getCardinalBPPARAM()

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

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

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