## ----style, echo=FALSE, results='asis'------------------------------------- BiocStyle::markdown() ## ----setup, echo=FALSE, message=FALSE--------------------------------------------------------------------------------- library(Cardinal) options(Cardinal.verbose=FALSE) options(Cardinal.progress=FALSE) options(width=120) 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, mean(.), .by="feature") # # # calculate tic # summarize(mse, tic=sum(.), .by="pixel") # # # calculate mean spectrum of circle region # mse %>% # select(region == "circle") %>% # summarize(mean(.)) ## ----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="median") %>% select(x==10, y==10, run=="run0") %>% process(plot=TRUE, par=list(main="Median interpolation", layout=c(2,1))) mse %>% reduceBaseline(method="locmin") %>% select(x==10, y==10, run=="run0") %>% process(plot=TRUE, par=list(main="Local minima")) ## ----reduceBaseline--------------------------------------------------------------------------------------------------- mse_pre <- reduceBaseline(mse_pre, method="median") ## ----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="simple") %>% process(plot=TRUE, par=list(main="Simple SD noise", layout=c(3,1))) mse_pre %>% select(x==10, y==10, run=="run0") %>% process() %>% peakPick(method="adaptive") %>% process(plot=TRUE, par=list(main="Adaptive SD noise")) mse_pre %>% select(x==10, y==10, run=="run0") %>% process() %>% peakPick(method="mad") %>% process(plot=TRUE, par=list(main="Adaptive MAD noise")) ## ----peakPick--------------------------------------------------------------------------------------------------------- mse_peaks <- peakPick(mse_pre, method="simple") ## ----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()