inst/doc/Cardinal-2-stats.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")

## ----eda-data, fig.height=3, fig.width=3, fig.align='center'------------------
set.seed(2020)
mse <- simulateImage(preset=2, npeaks=10, dim=c(20,20), sdnoise=0.5,
					peakheight=c(2,4), representation="centroid")

design <- makeFactor(circle=mse$circle, square=mse$square,
						bg=!(mse$circle | mse$square))

image(mse, design ~ x * y, key=TRUE)

## ----eda-image, fig.height=2, fig.width=7-------------------------------------
image(mse, feature=c(1,4,7), layout=c(1,3))

## ----pca----------------------------------------------------------------------
pca <- PCA(mse, ncomp=3)

summary(pca)

## ----pca-image, fig.height=2, fig.width=7-------------------------------------
image(pca, values="scores", superpose=FALSE, layout=c(1,3))

## ----pca-loadings, fig.height=2, fig.width=7----------------------------------
plot(pca, values="loadings", superpose=FALSE, layout=c(1,3), lwd=2)

## ----pca-scores, fig.height=2.5, fig.width=4, fig.align='center'--------------
pca_scores <- DataFrame(resultData(pca, 1, "scores"))

plot(pca_scores, PC1 ~ PC2, groups=design, pch=20)

## ----colocalized--------------------------------------------------------------
coloc <- colocalized(mse, mz=1023)
coloc

## ----colocalized-images, fig.height=2, fig.width=7----------------------------
image(mse, mz=coloc$mz[1:3], layout=c(1,3))

## ----ssc-clustering-----------------------------------------------------------
ssc <- spatialShrunkenCentroids(mse, r=1, k=5, s=c(0,3,6,9))

summary(ssc)

## ----ssc-image, fig.height=3, fig.width=3-------------------------------------
image(ssc, model=list(s=9), values="probability")

## ----ssc-statistic, fig.height=2.5, fig.width=4, fig.align='center'-----------
plot(ssc, model=list(s=9), values="statistic", lwd=2)

## ----ssc-top------------------------------------------------------------------
ssc_top <- topFeatures(ssc, model=list(s=9), class == 1)
ssc_top

## ----dgmm---------------------------------------------------------------------
dgmm <- spatialDGMM(mse, r=1, k=5, method="adaptive")

summary(dgmm)

## ----dgmm-image, fig.height=2, fig.width=7------------------------------------
image(dgmm, model=list(feature=c(1,4,7)), layout=c(1,3))

## ----dgmm-plot, fig.height=2, fig.width=7-------------------------------------
plot(dgmm, model=list(feature=c(1,4,7)), layout=c(1,3))

## ----dgmm-colocalized---------------------------------------------------------
coloc2 <- colocalized(dgmm, mse$square)
subset(coloc2, select=c(-r, -k, -group))

## ----classification-data, fig.height=2, fig.width=7---------------------------
set.seed(2020)
mse2 <- simulateImage(preset=7, npeaks=10, dim=c(10,10), sdnoise=0.5,
					nruns=3, peakdiff=2, representation="centroid")

class <- makeFactor(A=mse2$circleA, B=mse2$circleB)

image(mse2, class ~ x * y, key=TRUE, layout=c(1,3))

## ----classification-images, fig.height=2, fig.width=7-------------------------
image(mse2, feature=1, layout=c(1,3))

## ----pls-cv-------------------------------------------------------------------
cv_pls <- crossValidate(mse2, .y=class, .fun=PLS, ncomp=1:5, .fold=run(mse2))

summary(cv_pls)

## ----pls----------------------------------------------------------------------
pls <- PLS(mse2, y=class, ncomp=3)

summary(pls)

## ----pls-image, fig.height=2, fig.width=7-------------------------------------
image(pls, values="fitted", layout=c(1,3))

## ----pls-coefficients, fig.height=2.5, fig.width=4, fig.align='center'--------
plot(pls, values="coefficients", lwd=2)

## ----pls-scores, fig.height=2.5, fig.width=4, fig.align='center'--------------
pls_scores <- DataFrame(resultData(pls, 1, "scores"))

plot(pls_scores, C1 ~ C2, groups=class, pch=20)

## ----ssc-cv-------------------------------------------------------------------
cv_ssc <- crossValidate(mse2, .y=class,
						.fun=spatialShrunkenCentroids,
						r=1, s=c(0,3,6,9), .fold=run(mse2))

summary(cv_ssc)

## ----ssc-classification-------------------------------------------------------
ssc2 <- spatialShrunkenCentroids(mse2, y=class, r=1, s=0)

summary(ssc2)

## ----ssc-image-2, fig.height=2, fig.width=7-----------------------------------
image(ssc2, layout=c(1,3))

## ----ssc-statistic-2, fig.height=2.5, fig.width=4, fig.align='center'---------
plot(ssc2, values="statistic", lwd=2)

## ----ssc-top-2----------------------------------------------------------------
topFeatures(ssc2, class=="B") %>% subset(select=c(-diff, -r, -k, -s))

## ----test-data, fig.height=5, fig.width=7-------------------------------------
set.seed(2020)
mse3 <- simulateImage(preset=4, npeaks=10, dim=c(10,10), sdnoise=0.3,
					nruns=3, peakdiff=1, representation="centroid")

trt <- makeFactor(A=mse3$circleA, B=mse3$circleB)

image(mse3, trt ~ x * y, key=TRUE, layout=c(2,3))

## ----test-image, fig.height=5, fig.width=7------------------------------------
image(mse3, feature=1, layout=c(2,3))

## ----test-diff----------------------------------------------------------------
featureData(mse3)

## ----test-mean-test-----------------------------------------------------------
mtest <- meansTest(mse3, ~ condition, groups=run(mse3))

summary(mtest)

## ----test-mean-plot, fig.height=3, fig.width=7--------------------------------
plot(mtest, layout=c(2,5), ylab="intensity")

## ----test-mean-top------------------------------------------------------------
topFeatures(mtest, p.adjust="fdr", AdjP < .05)

## ----test-segment-dgmm--------------------------------------------------------
dgmm2 <- spatialDGMM(mse3, r=1, k=5, groups=run(mse3))

## ----test-segment-test--------------------------------------------------------
stest <- segmentationTest(dgmm2, ~ condition, classControl="Ymax")

summary(stest)

## ----test-segment-plot, fig.height=3, fig.width=7-----------------------------
plot(stest, layout=c(2,5), ylab="intensity")

## ----test-segment-mapping-----------------------------------------------------
image(stest, model=list(feature=3), values="mapping")

## ----test-segment-top---------------------------------------------------------
topFeatures(stest, p.adjust="fdr", AdjP < .05) %>% subset(select=c(-diff, -r, -k))

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