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

## ----setup, echo=FALSE, message=FALSE--------------------------------------
library(CardinalWorkflows)
register(SerialParam())
options(Cardinal.verbose=FALSE)
options(Cardinal.progress=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 <- summarize(rcc, .stat="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 <- summarize(rcc_peaks, .stat="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()