## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  error = FALSE,
  collapse = TRUE,
  comment = "#>"
)

## ----get-data-----------------------------------------------------------------
library(metabom8)

## ----import, fig.width=7.2, fig.height=4--------------------------------------
## Example: import 1D spectra
exp_dir <- system.file("extdata", package = "metabom8")
exp_type <- list(pulprog="noesygppr1d")

nmr_data <- read1d_proc(exp_dir, exp_type, n_max = 100)

names(nmr_data)

## ----pars_table, echo=FALSE---------------------------------------------------
knitr::kable(
  data.frame(
    Prefix = c("a_", "p_"),
    Meaning = c("Acquisition parameters",
                "Processing parameters"),
    Examples = c("a_PULPROG, a_SFO1, a_RG, a_SW_h, a_TE",
                 "p_WDW, p_LB, p_SI, p_PHC*, p_BC*, p_NC_proc")
  ),
  align = "lll"
)

## ----read-in-raw, fig.show='hold'---------------------------------------------
# import 1D NMR data
nmr_raw <- read1d_raw(
  exp_dir,
  exp_type,
  apodisation = list(fun = "exponential", lb = 0.2),
  zerofil = 1L,
  mode = "absorption"
  )

names(nmr_raw)

## ----apod, eval=TRUE, fig.width=8---------------------------------------------
# selected FID apodisation functions
f_apod <- c("sine","cosine","exponential","sem")

# create toy fid
n <- 200; t <- seq(0,1,len=n)
fid <- exp(-5*t)*cos(20*pi*t)

# apply apodisation
A  <- sapply(f_apod, \(f) metabom8:::.fidApodisationFct(n, list(fun=f, lb=-0.2)))
Fp <- sweep(A, 1, fid, `*`)
S  <- apply(Fp, 2, \(x) Mod(fft(x))[1:(n/2)])

# compare graphically
cols <- 1:ncol(A)

par(mfrow=c(2,2), mar=c(3,3,2,1))
plot(t, fid, type="l", lwd=2, main="FID")
matplot(t, A, type="l", lwd=2, col=cols, main="Apodisation windows")
matplot(t, Fp, type="l", lwd=2, col=cols, main="Windowed FID")
matplot(S, type="l", lwd=2, col=cols, main="Spectra")

legend("topright", f_apod, col=cols, lty=1, bty="n", cex=.8)

## ----viz_covid, fig.width=7.2, fig.height=4-----------------------------------
data("covid", package = "metabom8")

## plot directly from the dataset object
plot_spec(covid, shift = c(0.5, 2))

## alternatively extract data components: X, ppm
X <- covid$X
ppm <- covid$ppm

dim(X)
head(ppm)

plot_spec(X[1:3, ], ppm, shift = c(0.5, 2)) # backend='plotly' by default
plot_spec(X[1:3, ], ppm, shift = c(0.5, 2), backend='base')
plot_spec(X[1:3, ], ppm, shift = c(0.5, 2), backend='ggplot2')

## ----preproc_pipe-------------------------------------------------------------
data("hiit_raw", package = "metabom8") # urine NMR (HIIT exercise experiment)

names(hiit_raw) # X, ppm, meta

# piped preprocessing 
hiit_proc <- hiit_raw |> 
  calibrate(type = "tsp") |>
  excise() |>
  correct_baseline(method='asls') |>
  align_spectra() |>
  pqn()

dim(hiit_proc$X)

## ----viz-raw-hiit, fig.width=7.2, fig.height=4--------------------------------
plot_spec(hiit_raw, shift = c(4,4.1))

## ----viz-proc-hiit, fig.width=7.2, fig.height=4-------------------------------
plot_spec(hiit_proc, shift = c(4,4.1))

## ----preproc_stepwise, fig.width=7.2, fig.height=4----------------------------
X <- hiit_raw$X
ppm <- hiit_raw$ppm

# perform TSP calibration
X_cal <- calibrate(X, ppm, type='tsp')

# excise chemical shift regions
regions <- list(
      upfield_noise   = c(min(ppm), 0.25),
      residual_water  = c(4.5, 5.2),
      downfield_noise = c(9.7, max(ppm))
    )
X_exi <- excise(X_cal, ppm, regions)
ppm_exi <- attr(X_exi, 'm8_axis')$ppm

# baseline correction
X_bli <- correct_baseline(X_exi)

# plot spectrum
plot_spec(X_bli[1,], ppm_exi, shift = c(0,10))

## ----preproc_provenance-------------------------------------------------------

## Inspect preprocessing steps
print_provenance(hiit_proc)

# Retrieve specific step
get_provenance(hiit_proc, step = 2)

## Retrieve the full provenance log
log <- get_provenance(hiit_proc)
# print(log) # very detailed

## ----preproc_provenance-indi--------------------------------------------------
# Retrieve parameter from a named step
get_provenance(hiit_proc, step = "calibrate", param = "target")

## ----provenance_note----------------------------------------------------------
# collect env varaibles
params <- list(
  agent    = paste0("snakemake/", Sys.getenv("SNAKEMAKE_VERSION", "v?")),
  runtime  = "docker",
  workflow = Sys.getenv("M8_WORKFLOW", "std_prof-urine"),
  run_id   = Sys.getenv("M8_RUN_ID", "m8-2605-001")
)
 
# create note with env parameters
hiit_proc <- add_note(hiit_proc, 
         'Excluding outlier sample XX - reason is ...', 
         params)

# check provenance out
get_provenance(hiit_proc , step = "user note")

## ----mva-data-----------------------------------------------------------------
Y <- covid$an$type
X <- covid$X

## ----model_setup--------------------------------------------------------------
# Define the scaling and centering strategy
uv = uv_scaling(center = TRUE)

## ----pca----------------------------------------------------------------------
pca_model <- pca(X, scaling = uv, ncomp=2)

# show / summary methods 
show(pca_model)
summary(pca_model)

## ----resampling---------------------------------------------------------------
mc_cv <- balanced_mc(k=15, split=2/3, type="DA")

## ----pls----------------------------------------------------------------------
pls_model <- pls(X, Y, validation_strategy=mc_cv, scaling = uv)

# model information
show(pls_model)

# model performance summary
summary(pls_model)

# model scores & loadings
Tp <- scores(pls_model)
Pp <- loadings(pls_model)

## ----opls---------------------------------------------------------------------
opls_model <- opls(X, Y, validation_strategy=mc_cv, scaling = uv)

# model information & performance
show(opls_model)

## ----opls-summary-------------------------------------------------------------
summary(opls_model)

## ----opls-t-p-----------------------------------------------------------------
# model scores, loadings & vip
Tp <- scores(opls_model)
To <- scores(opls_model, orth=TRUE)

Pp <- loadings(opls_model)
Po <- loadings(opls_model, orth=TRUE)

vip <- vip(opls_model)

## ----opls-dx, fig.width=7.2, fig.height=4-------------------------------------
# distance to the model in X space
dex <- dmodx(opls_model)
head(dex)

## ----citation, echo=FALSE, results='asis'-------------------------------------
citation("metabom8")

## ----session-info, echo=FALSE-------------------------------------------------
sessionInfo()