## ----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=200) ## ----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) ## ----------------------------------------------------------------------------- spectra(mse) <- as.matrix(spectra(mse)) imageData(mse) ## ----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) ## ----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()