WiNN Tutorial: A Reproducible LC-MS Example

The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.

This vignette runs a single, fully reproducible WiNN workflow on simulated LC-MS data. Every object is created inside the document, so the example can be rerun unchanged during R CMD check. To keep the vignette lightweight, it uses only base R plus winn, avoids slow parameter-grid searches, and focuses on a single before/after comparison.

1. Simulate a small LC-MS dataset

We simulate 30 metabolites measured across 96 injections. Every eighth injection is a pooled QC sample, so QC values should be close after technical correction. The observed matrix includes dilution, drift, batch shifts, and random noise layered on top of a known ground-truth signal.

set.seed(42)

n_met <- 30L
n_samples <- 96L
n_batches <- 4L

batch <- rep(seq_len(n_batches), each = n_samples / n_batches)
run_order <- seq_len(n_samples)
qc_idx <- seq(4, n_samples, by = 8)
study_idx <- setdiff(seq_len(n_samples), qc_idx)

base_mean <- rnorm(n_met, mean = 7.8, sd = 0.35)
true_log <- base_mean + matrix(rnorm(n_met * n_samples, sd = 0.08), nrow = n_met)

signal_metabolites <- 1:6
study_signal <- as.numeric(scale(sin(study_idx / 9) + rnorm(length(study_idx), sd = 0.2)))
true_log[signal_metabolites, study_idx] <- true_log[signal_metabolites, study_idx] +
  0.30 * matrix(
    study_signal,
    nrow = length(signal_metabolites),
    ncol = length(study_idx),
    byrow = TRUE
  )

covariate_metabolites <- 7:10
continuous_covariate <- numeric(n_samples)
continuous_covariate[study_idx] <- as.numeric(scale(runif(length(study_idx))))
true_log[covariate_metabolites, ] <- true_log[covariate_metabolites, ] +
  0.20 * matrix(
    continuous_covariate,
    nrow = length(covariate_metabolites),
    ncol = n_samples,
    byrow = TRUE
  )

pooled_qc <- rowMeans(true_log[, study_idx, drop = FALSE])
true_log[, qc_idx] <- pooled_qc

make_batch_drift <- function(n, amplitude) {
  x <- seq(0, 1, length.out = n)
  amplitude * (0.8 * x - 0.6 * x^2 + 0.15 * sin(2 * pi * x))
}

batch_drift <- numeric(n_samples)
drift_amplitude <- c(0.22, 0.14, 0.18, 0.10)
for (b in seq_len(n_batches)) {
  idx <- which(batch == b)
  batch_drift[idx] <- make_batch_drift(length(idx), drift_amplitude[b])
}

metabolite_drift_scale <- matrix(
  0.05 + 0.25 * runif(n_met),
  nrow = n_met,
  ncol = n_samples
)
drift_log <- matrix(batch_drift, nrow = n_met, ncol = n_samples, byrow = TRUE) *
  metabolite_drift_scale

batch_shift <- matrix(0, nrow = n_met, ncol = n_samples)
batch_centers <- c(-0.18, 0.05, 0.12, -0.08)
for (b in seq_len(n_batches)) {
  idx <- which(batch == b)
  batch_shift[, idx] <- rnorm(n_met, mean = batch_centers[b], sd = 0.04)
}

dilution_factor <- exp(rnorm(n_samples, sd = 0.06))
noise_log <- matrix(rnorm(n_met * n_samples, sd = 0.08), nrow = n_met)

observed_log <- true_log +
  matrix(log(dilution_factor), nrow = n_met, ncol = n_samples, byrow = TRUE) +
  drift_log +
  batch_shift +
  noise_log

true_intensity <- pmax(expm1(true_log), 0)
observed_intensity <- pmax(expm1(observed_log), 0)

2. Run WiNN

For a vignette, fixed settings are the most practical choice because they are fast and deterministic. The auto-tuning mode is still available for full analyses when you want QC-guided parameter selection.

corrected_intensity <- winn(
  observed_intensity,
  batch = batch,
  run_order = run_order,
  control_samples = qc_idx,
  parameters = "fixed",
  fdr_threshold = 0.05,
  median_adjustment = "shrink",
  remove_batch_effects = "anova",
  lag = NULL,
  scale_by_batch = FALSE
)

3. Summarize agreement with the known truth

We use three metrics:

Mean per-metabolite Pearson correlation on the log1p scale
Mean per-metabolite ICC(A,1) on the log1p scale
Mean QC coefficient of variation (QC-CV) on the intensity scale

Higher Pearson and ICC indicate stronger agreement with the ground-truth profiles. Lower QC-CV indicates tighter pooled controls.

mean_qc_cv <- function(x, qc_samples) {
  qc_values <- x[, qc_samples, drop = FALSE]
  qc_cvs <- apply(qc_values, 1, function(v) {
    mu <- mean(v, na.rm = TRUE)
    if (!is.finite(mu) || abs(mu) < .Machine$double.eps) {
      return(NA_real_)
    }
    sd(v, na.rm = TRUE) / abs(mu)
  })
  mean(qc_cvs, na.rm = TRUE)
}

icc_a1 <- function(reference, candidate) {
  keep <- is.finite(reference) & is.finite(candidate)
  reference <- reference[keep]
  candidate <- candidate[keep]
  n <- length(reference)
  k <- 2L

  if (n < 2L) {
    return(NA_real_)
  }

  ratings <- cbind(reference, candidate)
  grand_mean <- mean(ratings)
  row_means <- rowMeans(ratings)
  col_means <- colMeans(ratings)

  ss_rows <- k * sum((row_means - grand_mean)^2)
  ss_cols <- n * sum((col_means - grand_mean)^2)
  ss_total <- sum((ratings - grand_mean)^2)
  ss_error <- ss_total - ss_rows - ss_cols

  ms_rows <- ss_rows / (n - 1)
  ms_cols <- ss_cols / (k - 1)
  ms_error <- ss_error / ((n - 1) * (k - 1))

  denom <- ms_rows + (k - 1) * ms_error + (k * (ms_cols - ms_error) / n)
  if (!is.finite(denom) || denom <= 0) {
    return(NA_real_)
  }

  (ms_rows - ms_error) / denom
}

calc_metrics <- function(candidate, truth, qc_samples) {
  candidate_log <- log1p(candidate)
  truth_log <- log1p(truth)

  pearson_vals <- vapply(seq_len(nrow(candidate)), function(i) {
    cor(candidate_log[i, ], truth_log[i, ], use = "pairwise.complete.obs")
  }, numeric(1))

  icc_vals <- vapply(seq_len(nrow(candidate)), function(i) {
    icc_a1(truth_log[i, ], candidate_log[i, ])
  }, numeric(1))

  data.frame(
    mean_pearson = mean(pearson_vals, na.rm = TRUE),
    mean_icc = mean(icc_vals, na.rm = TRUE),
    mean_qc_cv = mean_qc_cv(candidate, qc_samples)
  )
}

summary_metrics <- rbind(
  raw = calc_metrics(observed_intensity, true_intensity, qc_idx),
  winn = calc_metrics(corrected_intensity, true_intensity, qc_idx)
)

knitr::kable(
  round(summary_metrics, 4),
  caption = "Before/after metrics for the simulated dataset."
)

Before/after metrics for the simulated dataset.
	mean_pearson	mean_icc	mean_qc_cv
raw	0.5778	0.4843	0.1810
winn	0.6249	0.5577	0.1237

On this simulated dataset, WiNN should increase agreement with the known truth while reducing pooled-QC variability.

4. Visualize a few metabolites

The plot below shows four representative metabolites. Grey traces are the observed measurements, dashed black traces are the true underlying signal, and blue traces are the WiNN-corrected values. Orange points mark pooled QC injections.

met_ids <- c(1, 6, 9, 14)
panel_breaks <- (n_samples / n_batches) * seq_len(n_batches - 1) + 0.5

op <- par(mfrow = c(2, 2), mar = c(3.5, 3.5, 2.5, 1))

for (met in met_ids) {
  y_raw <- log1p(observed_intensity[met, ])
  y_true <- log1p(true_intensity[met, ])
  y_winn <- log1p(corrected_intensity[met, ])
  y_lim <- range(c(y_raw, y_true, y_winn), finite = TRUE)

  plot(
    run_order,
    y_raw,
    type = "l",
    col = "grey55",
    lwd = 1,
    ylim = y_lim,
    xlab = "Run order",
    ylab = "log1p(Intensity)",
    main = paste("Metabolite", met)
  )
  abline(v = panel_breaks, lty = 3, col = "grey80")
  lines(run_order, y_true, col = "black", lty = 2, lwd = 1.2)
  lines(run_order, y_winn, col = "#1b6ca8", lwd = 1.4)
  points(qc_idx, y_winn[qc_idx], pch = 16, cex = 0.6, col = "#d95f02")
}


par(op)

5. Takeaways

The vignette is self-contained: it generates its own data and reruns without external files or internet access.
A fixed-parameter WiNN run is fast enough for documentation while still showing the full correction workflow.
Pearson, ICC, and QC-CV provide complementary checks: signal recovery, agreement, and QC tightness.

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.