---
title: "Getting Started"
author: "Torben Kimhofer"
date: "`r Sys.Date()`"
output: 
  BiocStyle::html_document:
    toc: true
    number_sections: true
vignette: >
  %\VignetteIndexEntry{Getting Started}
  %\VignetteEncoding{UTF-8}
  %\VignetteEngine{knitr::rmarkdown}
editor_options: 
  markdown: 
    wrap: 72
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  error = FALSE,
  collapse = TRUE,
  comment = "#>"
)
set.seed(123)
```

# Introduction

This vignette demonstrates data import and preprocessing workflow for NMR-based 
metabolomics with **metabom8**. 

```{r get-data}
library(metabom8)
```

# Importing spectral data

In routine high-throughput NMR workflows, raw time-domain data (the
free-induction decay) are typically processed directly on the
spectrometer using the vendor software TopSpin. Downstream metabolomics
analysis therefore commonly starts from frequency-domain spectra. Consequently, 
this vignette focuses on importing spectra generated by TopSpin.

For completeness and teaching purposes, the section
([Starting point: raw FID](#import-fid)) briefly illustrates importing and 
transforming time domain data to obtain spectra.


## TopSpin-processed spectra {#import-spectra}

The function `read1d_proc()` imports TopSpin-processed 1D NMR
spectra from the `pdata` subdirectory of Bruker experiment folders
(typically `pdata/1`). It reads the processed absorption-mode spectrum
(e.g. `1r`) together with associated acquisition (`acqus`) and processing
(`procs`) parameters.

At minimum two arguments are required:

-   `path`: path to the directory that encloses NMR experiments
-   `exp_type`: list of spectrometer parameters to filter desired
    experiment types based on acquisition and processing metadata.
    
Further details on supported filtering options are provided in `?read1d_proc`.

```{r 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)
```

Both import functions (`read1d_proc()` and `read1d_raw()`) return a named
list containing three core components used throughout this vignette:
the spectral matrix (`X`), the chemical-shift axis (`ppm`), and associated
metadata (`meta`). 

For interactive exploration, setting `to_global = TRUE` assigns the
components `X`, `ppm`, and `meta` directly to the global environment
(`.GlobalEnv`). Existing objects with these names will be overwritten.

The row names of `X` and `meta` correspond to experiment directories and can be
used to join external sample annotations. The `meta` data frame
contains full TopSpin acquisition and processing parameters for each
spectrum.

Example acquisition and processing parameters as they appear in `meta`:

```{r 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"
)
```

See `?read1d_proc()` for further details.

## Raw FID data {#import-fid}

For completeness, `metabom8` also provides functionality to import and
process raw FIDs. This allows users to reproduce basic spectral processing
steps within R or to explore how different processing parameters influence
the resulting spectra.

Processing an FID converts the time-domain signal into a frequency-domain
spectrum. Typical processing steps include digital-filter correction 
(group-delay), apodisation (windowing), zero-filling, Fourier 
transformation, phase correction.

The function `read1d_raw()` performs these operations:

```{r 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)
```

The following examples illustrate several commonly used apodisation
profiles implemented internally in `metabom8`.

```{r 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)
```

## Visual inspection

Spectra can be visualised with the function `plot_spec()`. The minimum
set of arguments include the spectral data `X` (array or matrix) and the
chemical shift vector `ppm`.

Three plotting backends are available: `"plotly"` (default), `"base"`, and
`"ggplot2"`. The `"plotly"` backend is set up to use WebGL, enabling
smooth interactive graphics even when the number of data spectra is
large. The `"ggplot2"` backend is not iteractive and renders much slower. 
However, the returned object integrates with the extensive `"ggplot2"` ecosystem 
and can be further customised to generate publication-quality figures.

See `?plot_spec` for additional details.

```{r 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')
```

# Preprocessing

`metabom8` provides a modular preprocessing API for 1D NMR spectra. Each
preprocessing function is composable and designed to support both
interactive exploration and reproducible pipelines.

For a list of preprocessing functions call `list_preprocessing()` from
the R console.

## Pipeline-style usage

When operating on a `metabom8`-style dataset object, preprocessing functions return
the same structured object: a named list containing 
the elements `X`, `ppm`, and `meta`. This allows steps to be chained using the
base R pipe operator (`|>`):

```{r 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)
```

This API design enables clear, declarative workflows suitable for production 
analyses.



```{r viz-raw-hiit, fig.width=7.2, fig.height=4}
plot_spec(hiit_raw, shift = c(4,4.1))
```

```{r viz-proc-hiit, fig.width=7.2, fig.height=4}
plot_spec(hiit_proc, shift = c(4,4.1))
```

## Stepwise execution

Each operation can also be applied independently to a spectral matrix
and ppm vector. This explicit alternative is useful for pipeline development,
parameter tuning and visual inspection.

```{r 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))
```

## Provenance

Data provenance refers to the documented record of how data were generated,
processed, and transformed. Recording this information ensures that analyses
remain reproducible and that individual preprocessing steps can be inspected,
verified, or repeated at a later stage. 

In `metabom8`, particular emphasis is placed on transparent preprocessing, as 
these decisions can substantially influence downstream statistical analyses and 
biological interpretation. Each preprocessing operation automatically records 
its parameters and execution context as structured metadata attached to the 
spectral matrix `X`.

The resulting dataset therefore remains self-describing: the complete processing 
history is linked to the data matrix and can be inspected programmatically at 
any time.

```{r 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
```

Each entry stores: 

* sequential execution index
* preprocessing step name 
* parameter settings
* optional notes
* timestamp and versioning info

Individual parameter values can be retrieved as well, allowing the full transformation 
history of a dataset to be inspected programmatically.

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

Custom notes can also be included in the provenance log. This is particularly userful
when data are processed in automated fashion, for example when using workflow managers:

```{r 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")
```


> **Note on provenance**  
`metabom8` records preprocessing history as structured attributes on
standard R matrices or named lists. This design is lightweigth and prioritises
interoperability and flexibility over strict container classes
(e.g., `SummarizedExperiment`). Provenance is preserved when objects are saved/loaded
based on standard R serialisation (`save()`/`saveRDS()`, `load()`), but may be lost 
if the object is modified by functions outside the `metabom8`
API that drop or overwrite attributes (for example when objects are
coerced, reconstructed, or exported to external formats).


# Modelling

The section demonstrates data modelling using PCA (unsupervised) and PLS/OPLS (supervised).

All modelling calls return an instance of class `m8_model`, which can be used 
for further processing and diagnostics.


```{r mva-data}
Y <- covid$an$type
X <- covid$X
```

To prevent the natural dynamic range of metabolite concentrations from 
dominating statistical modelling, data are scaled prior to analysis. In 
**metabom8**, scaling and centering are passed explicitly to the model function:

```{r model_setup}
# Define the scaling and centering strategy
uv = uv_scaling(center = TRUE)
```

See `?uv_scaling` for further details and additional scaling strategies.



## Principal Components Analysis (PCA)

Unsupervised modelling with PCA is done by an iterative procedure
(NIPALS) that calculates components iteratively, one at a time. Therefore, the
number of desired components is provided in the function call.

Here we compute two principal components using unit-variance scaling:

```{r pca}
pca_model <- pca(X, scaling = uv, ncomp=2)

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

## Partial-Least Squares (PLS)
Partial-Least-Squares (PLS) as supervised modelling technique utilises 
a statistical re-sampling techniques to determine the optimal number of components 
and to safeguard against overfitting.

In this vignette, we employ balanced Monte Carlo cross-validation (`balanced_mc`). 
This process repeatedly splits the data into training and testing sets to ensure the 
model predictive power generalizes well to unseen data.

```{r resampling}
mc_cv <- balanced_mc(k=15, split=2/3, type="DA")
```

See `?balanced_mc` for further details and additional statistical validation strategies.

```{r 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)
```

## Orthogonal-Partial-Least Squares (OPLS)

OPLS can be though of as PLS with an integrated orthogonality filter, which has 
been reported to improve model interpretability. While standard 
PLS mixes predictive and non-predictive variation, OPLS partitions the variance 
into two parts: variation correlated to the response ($Y$) and variation "orthogonal" 
(unrelated) to it.

Similar to PLS, the OPLS algorithm evaluates model performance 
using statistical re-sampling techniques to ensure the results are robust and 
not driven by noise.

```{r opls}
opls_model <- opls(X, Y, validation_strategy=mc_cv, scaling = uv)

# model information & performance
show(opls_model)
```
The OPLS algorithm fitted one predictive and one orthogonal component.
A second orthogonal component was evaluated but the improvement in
cross-validated class prediction (`AUCv`)based on the left-out sets was
negligible (`stop reason`). Therefore, a third component was not fitted to the full
spectral matrix.

```{r opls-summary}
summary(opls_model)
```

For detailed information on component stopping criteria used in `opls()` and`pls()` 
see `?.evalFit()`.

Methods implemented for both, `pls()` and `opls()`, models include `scores()`, 
`loadings()` and `fitted()` (Y-fitted values from statistical resampling).
For OPLS models, these former two functions accept argument `orth = TRUE`
to extract orthogonal component information.

Further OPLS model interpretation can be performed using Variable Importance in 
Projection (`vip()`), which identifies the most influential variables in the model.

```{r 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)
```


To ensure the model is statistically robust, `metabom8` provides several 
diagnostic tools:

* `dmodx()`: Calculates the "Distance to the Model X" to identify spectral outliers
* `opls_perm()`:  Executes Y-permutations to validate that the observed $Q^2$ is 
significantly higher than what would be expected by chance.
* `cv_anova()`: Performs a Cross-Validated ANOVA to test the significance of the 
model's predictive ability (only for regression)


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

This vignette demonstrated the core data import, preprocessing, and modelling workflow 
implemented in `metabom8`. For further details on individual functions, see the 
corresponding help pages.

# Citation

```{r citation, echo=FALSE, results='asis'}
citation("metabom8")
```

# Session Info

```{r session-info, echo=FALSE} 
sessionInfo()
```
