## ----style, echo=FALSE, results='asis'---------------------------------------- BiocStyle::markdown() ## ----setup, echo=FALSE, message=FALSE----------------------------------------- library(Cardinal) ## ----install, eval=FALSE------------------------------------------------------ # install.packages("BiocManager") # BiocManager::install("Cardinal") ## ----library, eval=FALSE------------------------------------------------------ # library(Cardinal) ## ----read-imzML-continuous---------------------------------------------------- path_continuous <- CardinalIO::exampleImzMLFile("continuous") path_continuous mse_tiny <- readMSIData(path_continuous) mse_tiny ## ----read-imzML-processed----------------------------------------------------- path_processed <- CardinalIO::exampleImzMLFile("processed") path_processed msa_tiny <- readMSIData(path_processed) msa_tiny ## ----show-MSImagingArrays----------------------------------------------------- msa_tiny ## ----subset-MSImagingArrays--------------------------------------------------- msa_tiny[1:3] ## ----spectraData-MSImagingArrays---------------------------------------------- spectraData(msa_tiny) ## ----spectra-accessor--------------------------------------------------------- spectra(msa_tiny, "mz") spectra(msa_tiny, "intensity") ## ----intensity-accessor------------------------------------------------------- mz(msa_tiny) intensity(msa_tiny) ## ----pixelData-MSImagingArrays------------------------------------------------ pixelData(msa_tiny) pData(msa_tiny) ## ----coord-accessor----------------------------------------------------------- coord(msa_tiny) ## ----run-accessor------------------------------------------------------------- runNames(msa_tiny) head(run(msa_tiny)) ## ----show-MSImagingExperiment------------------------------------------------- mse_tiny ## ----subset-MSImagingExperiment----------------------------------------------- mse_tiny[1:500, 1:3] ## ----spectraData-MSImagingExperiment------------------------------------------ spectraData(mse_tiny) spectra(mse_tiny) ## ----pixelData-MSImagingExperiment-------------------------------------------- pixelData(mse_tiny) ## ----featureData-accessor----------------------------------------------------- featureData(mse_tiny) fData(mse_tiny) ## ----mz-accessor-------------------------------------------------------------- head(mz(mse_tiny)) ## ----constructor-------------------------------------------------------------- set.seed(2020, kind="L'Ecuyer-CMRG") s <- simulateSpectra(n=9, npeaks=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(spectraData=s$intensity, featureData=fdata, pixelData=pdata) out ## ----simulate----------------------------------------------------------------- # Simulate an MSImagingExperiment set.seed(2020, kind="L'Ecuyer-CMRG") mse <- simulateImage(preset=6, dim=c(32,32), baseline=0.5) mse # Create a version as MSImagingArrays msa <- convertMSImagingExperiment2Arrays(mse) msa ## ----plot-i, fig.height=3, fig.width=9---------------------------------------- plot(msa, i=c(496, 1520)) ## ----plot-coord, fig.height=3, fig.width=9------------------------------------ plot(msa, coord=list(x=16, y=16)) ## ----plot-superpose, fig.height=3, fig.width=9-------------------------------- plot(msa, i=c(496, 1520), xlim=c(1000, 1250), superpose=TRUE) ## ----image-i, fig.height=4, fig.width=9--------------------------------------- image(mse, i=1938) ## ----image-mz, fig.height=4, fig.width=9-------------------------------------- image(mse, mz=1003.3) ## ----image-plusminus, fig.height=4, fig.width=9------------------------------- image(mse, mz=1003.3, tolerance=0.5, units="mz") ## ----image-run, fig.height=4, fig.width=5------------------------------------- image(mse, mz=1003.3, run="runA1") ## ----image-subset, fig.height=4, fig.width=9---------------------------------- image(mse, mz=1003.3, subset=mse$circleA | mse$circleB) ## ----image-smooth, fig.height=4, fig.width=9---------------------------------- image(mse, mz=1003.3, smooth="gaussian") ## ----image-contrast, fig.height=4, fig.width=9-------------------------------- image(mse, mz=1003.3, enhance="histogram") ## ----image-superpose, fig.height=4, fig.width=9------------------------------- image(mse, mz=c(1003.3, 1663.6), superpose=TRUE, enhance="adaptive", scale=TRUE) ## ----select-ROI, eval=FALSE--------------------------------------------------- # sampleA <- selectROI(mse, mz=1003.3, subset=run(mse) == "run0") # sampleB <- selectROI(mse, mz=1003.3, subset=run(mse) == "run1") ## ----makeFactor, eval=FALSE--------------------------------------------------- # regions <- makeFactor(A=sampleA, B=sampleB) ## ----pdf, eval=FALSE---------------------------------------------------------- # pdffile <- tempfile(fileext=".pdf") # # pdf(file=pdffile, width=9, height=4) # image(mse, mz=1003.3) # dev.off() ## ----style-dark, fig.height=4, fig.width=9------------------------------------ image(mse, mz=1003.3, style="dark") ## ----print-------------------------------------------------------------------- p <- image(mse, mz=1003.3) plot(p, smooth="guide") ## ----subset-1----------------------------------------------------------------- # subset first 5 mass spectra msa[1:5] ## ----subset-2----------------------------------------------------------------- # subset first 10 images and first 5 mass spectra mse[1:10, 1:5] ## ----features----------------------------------------------------------------- # get row index corresponding to m/z 1003.3 features(mse, mz=1003.3) # get row indices corresponding to m/z 1002 - 1004 features(mse, 1002 < mz & mz < 1004) ## ----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 "runA1" pixels(mse, x <= 3, y <= 3, run == "runA1") ## ----subset-3----------------------------------------------------------------- fid <- features(mse, 1002 < mz, mz < 1004) pid <- pixels(mse, x <= 3, y <= 3, run == "runA1") mse[fid, pid] ## ----subset-method-1---------------------------------------------------------- # subset MSImagingArrays subset(msa, x <= 3 & x <= 3) ## ----subset-method-2---------------------------------------------------------- # subset MSImagingExperiment subset(mse, 1002 < mz & mz < 1004, x <= 3 & x <= 3) ## ----subsetFeatures----------------------------------------------------------- # subset features subsetFeatures(mse, 1002 < mz, mz < 1004) # subset pixels subsetPixels(mse, x <= 3, y <= 3) ## ----slice-------------------------------------------------------------------- # slice image for first mass feature a <- sliceImage(mse, 1) dim(a) ## ----slice-mz----------------------------------------------------------------- # slice image for m/z 1003.3 a2 <- sliceImage(mse, mz=1003.3, drop=FALSE) dim(a2) ## ----slice-image, fig.height=4, fig.width=9----------------------------------- par(mfcol=c(1,2), new=FALSE) image(a2[,,1,1], asp=1) image(a2[,,2,1], asp=1) ## ----cbind-divide------------------------------------------------------------- # divide dataset into separate objects for each run mse_run0 <- mse[,run(mse) == "runA1"] mse_run1 <- mse[,run(mse) == "runB1"] 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(A=mse$circleA, B=mse$circleB, other=mse$square1 | mse$square2) pData(mse) ## ----fData-set---------------------------------------------------------------- fData(mse)$region <- makeFactor( circle=mz(mse) > 1000 & mz(mse) < 1250, square=mz(mse) < 1000 | mz(mse) > 1250) fData(mse) ## ----spectra-set-------------------------------------------------------------- # create a new spectra matrix of log-intensities spectra(mse, "log2intensity") <- log2(spectra(mse) + 1) spectraData(mse) # examine the new spectra matrix spectra(mse, "log2intensity")[1:5, 1:5] ## ----centroided-get----------------------------------------------------------- centroided(mse) ## ----centroided-set, eval=FALSE----------------------------------------------- # centroided(mse) <- FALSE ## ----summarize-features, fig.height=3, fig.width=9---------------------------- # calculate mean spectrum mse <- summarizeFeatures(mse, stat="mean") # mean spectrum stored in featureData fData(mse) # plot mean spectrum plot(mse, "mean", xlab="m/z", ylab="Intensity") ## ----summarize-pixels, fig.height=4, fig.width=9------------------------------ # calculate TIC mse <- summarizePixels(mse, stat=c(TIC="sum")) # TIC stored in pixelData pData(mse) # plot TIC image(mse, "TIC", col=matter::cpal("Cividis")) ## ----summarize-features-groups, fig.height=3, fig.width=9--------------------- # calculate mean spectrum mse <- summarizeFeatures(mse, stat="mean", groups=mse$region) # mean spectrum stored in featureData fData(mse) # plot mean spectrum plot(mse, c("A.mean", "B.mean"), xlab="m/z", ylab="Intensity") ## ----summarize-pixels-groups, fig.height=8, fig.width=9----------------------- # calculate mean spectrum mse <- summarizePixels(mse, stat="sum", groups=fData(mse)$region) # mean spectrum stored in featureData pData(mse) # plot mean spectrum image(mse, c("circle.sum", "square.sum"), scale=TRUE) ## ----matter------------------------------------------------------------------- # spectra are stored as an out-of-memory matrix spectra(mse_tiny) spectraData(mse_tiny) # 'intensity' array is 'matter_mat' object ## ----------------------------------------------------------------------------- # coerce spectra to an in-memory matrix spectra(mse_tiny) <- as.matrix(spectra(mse_tiny)) spectraData(mse_tiny) # 'intensity' array is 'matrix' object ## ----coerce-1----------------------------------------------------------------- msa # coerce to MSImagingExperiment as(msa, "MSImagingExperiment") ## ----coerce-2----------------------------------------------------------------- mse # coerce to MSImagingArrays as(mse, "MSImagingArrays") ## ----normalize---------------------------------------------------------------- msa_pre <- normalize(msa, method="tic") ## ----smoothSignal-plot, fig.height=7, fig.width=9----------------------------- p1 <- smooth(msa, method="gaussian") |> plot(coord=list(x=16, y=16), xlim=c(1150, 1450), linewidth=2) p2 <- smooth(msa, method="sgolay") |> plot(coord=list(x=16, y=16), xlim=c(1150, 1450), linewidth=2) matter::as_facets(list(p1, p2), nrow=2, labels=c("Gaussian smoothing", "Savitsky-Golay smoothing")) ## ----smoothSignal------------------------------------------------------------- msa_pre <- smooth(msa_pre, method="gaussian") ## ----reduceBaseline-plot, fig.height=5, fig.width=9--------------------------- p1 <- reduceBaseline(mse, method="locmin") |> plot(coord=list(x=16, y=16), linewidth=2) p2 <- reduceBaseline(mse, method="median") |> plot(coord=list(x=16, y=16), linewidth=2) matter::as_facets(list(p1, p2), nrow=2, labels=c("Local minima interpolation", "Running medians")) ## ----reduceBaseline----------------------------------------------------------- msa_pre <- reduceBaseline(msa_pre, method="locmin") ## ----unaligned-spectra, fig.height=3, fig.width=9----------------------------- set.seed(2020, kind="L'Ecuyer-CMRG") mse_drift <- simulateImage(preset=1, npeaks=10, from=500, to=600, sdmz=750, units="ppm") plot(mse_drift, i=186:195, xlim=c(535, 570), superpose=TRUE, key=FALSE, linewidth=2) ## ----recalibrate, fig.height=3, fig.width=9----------------------------------- peaks_drift <- estimateReferencePeaks(mse_drift) mse_nodrift <- recalibrate(mse_drift, ref=peaks_drift, method="locmax", tolerance=1500, units="ppm") mse_nodrift <- process(mse_nodrift) plot(mse_nodrift, i=186:195, xlim=c(535, 570), superpose=TRUE, key=FALSE, linewidth=2) ## ----process-execute---------------------------------------------------------- msa_pre <- process(msa_pre) ## ----peakPick-plot, fig.height=7, fig.width=9--------------------------------- p1 <- peakPick(msa_pre, method="diff", SNR=3) |> plot(coord=list(x=16, y=16), linewidth=2) p2 <- peakPick(msa_pre, method="filter", SNR=3) |> plot(coord=list(x=16, y=16), linewidth=2) matter::as_facets(list(p1, p2), nrow=2, labels=c("Derivative-based SNR", "Dynamic filtering-based SNR")) ## ----peakPick----------------------------------------------------------------- msa_peaks <- peakPick(msa_pre, method="filter", SNR=3) ## ----peakAlign---------------------------------------------------------------- mse_peaks <- peakAlign(msa_peaks, tolerance=200, units="ppm") mse_peaks ## ----subsetFeatures-peaks----------------------------------------------------- fData(mse_peaks) # filter to peaks with frequencies > 0.1 mse_filt <- subsetFeatures(mse_peaks, freq > 0.1) fData(mse_filt) ## ----mean-peaks-1, fig.height=3, fig.width=9---------------------------------- mse_filt <- summarizeFeatures(mse_filt) plot(mse_filt, "mean", xlab="m/z", ylab="Intensity", linewidth=2, annPeaks=10) ## ----peakPick-ref------------------------------------------------------------- msa_peaks2 <- peakPick(msa_pre, ref=mz(mse_filt), type="area", tolerance=600, units="ppm") mse_peaks2 <- process(msa_peaks2) ## ----mean-peaks-2, fig.height=3, fig.width=9---------------------------------- mse_peaks2 <- as(mse_peaks2, "MSImagingExperiment") mse_peaks2 <- summarizeFeatures(mse_peaks2) plot(mse_peaks2, "mean", xlab="m/z", ylab="Intensity", linewidth=2, annPeaks=10) ## ----peakProcess-------------------------------------------------------------- mse_peaks3 <- peakProcess(msa_pre, method="diff", SNR=6, sampleSize=0.3, filterFreq=0.02) mse_peaks3 ## ----bin---------------------------------------------------------------------- mse_binned <- bin(msa, resolution=1, units="mz") mse_binned ## ----process-workflow, fig.height=5, fig.width=9------------------------------ mse_queue <- msa |> normalize() |> smooth() |> reduceBaseline() |> peakPick(SNR=6) # preview processing plot(mse_queue, coord=list(x=16, y=16), linewidth=2) # apply processing and align peaks mse_proc <- peakAlign(mse_queue) mse_proc <- subsetFeatures(mse_proc, freq > 0.1) mse_proc ## ----write-imzML-------------------------------------------------------------- imzfile <- tempfile(fileext=".imzML") writeMSIData(mse_proc, file=imzfile) list.files(imzfile) ## ----read-imzML--------------------------------------------------------------- mse_re <- readMSIData(imzfile) mse_re ## ----set-BPPARAM, eval=FALSE-------------------------------------------------- # # run in parallel, rather than serially # mse_mean <- summarizeFeatures(mse, BPPARAM=MulticoreParam()) ## ----registered--------------------------------------------------------------- BiocParallel::registered() ## ----getCardinalBPPARAM------------------------------------------------------- getCardinalBPPARAM() ## ----setCardinalBPPARAM-snow-------------------------------------------------- # register a SNOW backend setCardinalBPPARAM(SnowParam(workers=2, progressbar=TRUE)) getCardinalBPPARAM() ## ----setCardinalBPPARAM-null-------------------------------------------------- # reset backend setCardinalBPPARAM(NULL) ## ----RNGkind, eval=FALSE------------------------------------------------------ # set.seed(1, kind="L'Ecuyer-CMRG") ## ----install-CardinalWorkflows, eval=FALSE------------------------------------ # BiocManager::install("CardinalWorkflows") ## ----library-CardinalWorkflows, eval=FALSE------------------------------------ # library(CardinalWorkflows) ## ----session-info------------------------------------------------------------- sessionInfo()