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

## ----setup, echo=FALSE, message=FALSE-----------------------------------------
library(CardinalWorkflows)
setCardinalBPPARAM(SerialParam())
setCardinalVerbose(FALSE)
RNGkind("L'Ecuyer-CMRG")

## ----library------------------------------------------------------------------
library(Cardinal)

## ----load-rcc-----------------------------------------------------------------
data(rcc, package="CardinalWorkflows")
rcc <- as(rcc, "MSImagingExperiment")

## ----show-rcc-----------------------------------------------------------------
rcc

## ----rcc-mean-----------------------------------------------------------------
rcc_mean <- summarizeFeatures(rcc, "mean")

## ----rcc-peak-process---------------------------------------------------------
rcc_ref <- rcc_mean %>%
  peakPick(SNR=3) %>%
  peakAlign(ref="mean",
            tolerance=0.5,
            units="mz") %>%
  peakFilter() %>%
  process()

## ----rcc-peak-bin-------------------------------------------------------------
rcc_peaks <- rcc %>%
  normalize(method="tic") %>%
  peakBin(ref=mz(rcc_ref),
          tolerance=0.5,
          units="mz") %>%
  process()

rcc_peaks

## ----rcc-split----------------------------------------------------------------
xcutoff<-c(35, 23, 28, 39, 29, 28, 44, 32)

rcc_peaks$rough_diagnosis <- factor("normal", level=c("cancer", "normal"))

for ( i in 1:nlevels(run(rcc_peaks)) ) {
  cur_run <- run(rcc_peaks) == runNames(rcc_peaks)[i]
  pData(rcc_peaks)$rough_diagnosis[cur_run & coord(rcc_peaks)$x < xcutoff[i]] <- "cancer"
}

rcc_peaks$groups <- interaction(run(rcc_peaks), rcc_peaks$rough_diagnosis)

## ----rcc-check, fig.height=10-------------------------------------------------
image(rcc_peaks, mz=810, groups=rough_diagnosis,
      contrast.enhance="histogram", layout=c(4,2))

## ----rcc-var------------------------------------------------------------------
rcc_var <- summarizeFeatures(rcc_peaks, "var", as="DataFrame")

plot(rcc_var, var ~ mz, main="variance")

## ----rcc-filter---------------------------------------------------------------
rcc_peaks2 <- rcc_peaks[rcc_var$var >= quantile(rcc_var$var, 0.8),]

## ----dgmm1--------------------------------------------------------------------
set.seed(1)
rcc_dgmm1 <- spatialDGMM(rcc_peaks2[16,], r=1, k=4, groups=1)

summary(rcc_dgmm1)

## ----dgmm1-plot, fig.height=10------------------------------------------------
image(rcc_dgmm1, layout=c(4,2))

## ----dgmm---------------------------------------------------------------------
set.seed(1)
rcc_dgmm <- spatialDGMM(rcc_peaks2, r=1, k=4, groups=rcc_peaks2$groups)

summary(rcc_dgmm)

## ----mtest--------------------------------------------------------------------
mtest <- meansTest(rcc_peaks2, ~ rough_diagnosis, groups=rcc_peaks2$groups)

summary(mtest)

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

## ----stest--------------------------------------------------------------------
stest <- segmentationTest(rcc_dgmm, ~ rough_diagnosis)

summary(stest)

## ----stest-top----------------------------------------------------------------
topFeatures(stest, p.adjust="fdr", AdjP < .1)

## ----stest-plot---------------------------------------------------------------
plot(stest, model=list(feature=16))

## ----top-image, fig.height=10-------------------------------------------------
image(rcc_peaks2, mz=885, layout=c(4,2),
      contrast.enhance="suppress",
      normalize.image="linear")

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