bio-metabolomics-statistical-analysis

What this skill does

## Version Compatibility

Reference examples tested with: ropls 1.34+, scipy 1.12+, statsmodels 0.14+, numpy 1.26+, pandas 2.1+, matplotlib 3.8+

Before using code patterns, verify installed versions match. If versions differ:
- Python: `pip show <package>` then `help(module.function)` to check signatures
- R: `packageVersion('<pkg>')` then `?function_name` to verify parameters

If code throws ImportError, AttributeError, or TypeError, introspect the installed
package and adapt the example to match the actual API rather than retrying.

# Metabolomics Statistical Analysis

**"Tell me which metabolites separate my groups"** -> run an honest univariate test with dependence-aware FDR AND a permutation-validated multivariate model, then reconcile the two.
- R: `ropls::opls()` (PCA/PLS-DA/OPLS-DA + permutation), `wilcox.test()`/`lm()`/`lme4::lmer()`, `p.adjust(method='BH')`
- Python: `scipy.stats.ttest_ind(equal_var=False)`/`mannwhitneyu`, `statsmodels` `multipletests(method='fdr_bh')`, `sklearn.cross_decomposition.PLSRegression` + `permutation_test_score`

## The Single Most Important Modern Insight -- A Score Plot Is a Hypothesis, Not a Result

In metabolomics the regime is p >> n (hundreds-to-thousands of features, tens of samples) with strongly correlated features. In that regime any binary labelling of n points in >= n-1 dimensions is linearly separable with probability 1, so PLS-DA and OPLS-DA produce a clean two-cluster score plot even for randomly assigned labels. A beautiful score plot is the generic output of the algorithm and carries essentially zero information. Only a cross-validated Q2 benchmarked against a permutation null distinguishes signal from geometry (Westerhuis 2008; Ruiz-Perez 2020). Three corollaries reorganize the whole skill: (1) R2 is no evidence (it can be driven to 1 by adding components); only permutation-validated Q2 licenses a claim. (2) Scaling is a hidden hypothesis -- variance-driven methods weight a feature by the variance it is allowed to contribute, so Pareto vs unit-variance hands back a different VIP list and a different biological story (van den Berg 2006). (3) Features are not independent (pathways, adducts, isotopologues), so naive BH-independence is violated and one real signal lights up its whole correlated cluster.

## Transformation and Scaling -- the Decision That Changes Conclusions

Transformation (nonlinear, per value: corrects heteroscedastic multiplicative MS noise and skew) and scaling (linear, per feature: sets relative weight) are distinct. Mean-centering is the universal first step. Choosing not to scale is the strongest prior of all -- it lets the most abundant metabolite drive PC1.

| Method | Per feature j | Effect | Use when |
|--------|---------------|--------|----------|
| Centering | subtract mean | offsets removed, variance unchanged | always (prerequisite for all below) |
| Auto / unit-variance / "standard" | center, / SD_j | every metabolite equal weight | a priori all metabolites equally important; classic default -- but inflates near-LOD noise |
| Pareto | center, / sqrt(SD_j) | between raw and UV | de facto metabolomics/NMR/OPLS-DA default; curbs dominance with less noise inflation than UV |
| Range | center, / (max-min) | abundance dependence removed | clean data, few outliers (outlier-sensitive) |
| Vast | UV x (mean/SD) | up-weights low-CV stable features | focus on robust/reproducible features with prior class info |
| Level | center, / mean_j | relative (% change) | response as relative change (mean noisy at low abundance) |
| Log / log10 | log(x) | multiplicative -> additive | concentrations spanning orders of magnitude (undefined at 0) |
| glog | linear near 0, log for large x | variance-stabilizing | data with zeros / near-LOD values; preferred over plain log (needs transition param) |
| Power (sqrt, cube-root) | x^(1/2) | mild stabilization | mild skew with zeros present |

van den Berg 2006: on real data autoscale and range recovered biologically meaningful loadings; Pareto is the pragmatic middle. Decision rule: transform first (if heteroscedastic -- usually yes for MS), then center, then scale; run at least Pareto AND UV, and if the top-VIP list or conclusion flips, the result is scaling-fragile and must be tempered.

## Decision Tree by Scenario

| Goal / situation | Do | Why |
|------------------|----|----|
| First look, QC, batch/outlier check | PCA on scaled data; color scores by batch/injection order; Hotelling T2 ellipse | Unsupervised -> cannot overfit the grouping; pooled-QC samples must cluster tightly in the center, else analytical variance dominates |
| Which single metabolites differ (2 groups) | Welch t-test (post-transform) or Mann-Whitney; BH FDR; report fold change + CI | Interpretable per-feature effect; FDR-controlled; effect size mandatory in p>>n |
| 2 groups, paired/pre-post | Paired t-test or Wilcoxon signed-rank | Discards within-subject pairing if analyzed unpaired -> underpowered |
| >2 groups | One-way ANOVA (+Tukey) or Kruskal-Wallis (+Dunn) | Match normality assumption |
| Longitudinal / repeated measures | Linear mixed model (random intercept/slope per subject) | Handles unbalanced timepoints, missingness, within-subject correlation |
| Covariate adjustment (age/sex/BMI/batch) | Per-feature linear model `y ~ group + covars` | Human metabolome is dominated by age/sex/BMI -- unadjusted, they masquerade as case/control signal (Thevenot 2015) |
| A discriminant / predictive model | PLS-DA (`orthoI=0`) or OPLS-DA (`predI=1, orthoI=NA`) + permutation + double CV | Supervised; demands full validation (see checklist) |
| Built-in feature selection | sparse PLS-DA `splsda()` (mixOmics) with `tune.splsda` | Selection must be inside CV -> hand off to machine-learning/biomarker-discovery |
| Rank discriminant features | VIP from a permutation-validated model only; corroborate with univariate FDR | VIP > 1 is a heuristic, not a test (see failure modes) |
| Confirm a biomarker | Independent validation cohort | Internal CV does not correct overfitting/forking-paths; discovery performance overestimates external |

## Scaling + PCA

**Goal:** Get the honest unsupervised first look that cannot chase the labels, with QC as the primary data-quality readout.

**Approach:** Transform if heteroscedastic, then PCA with an explicit scaling; inspect QC clustering, batch coloring, and Hotelling T2.

```r
library(ropls)
# scaleC default is "standard" (unit-variance/autoscale), NOT Pareto -- set explicitly
pca <- opls(t(feature_matrix), scaleC = 'pareto', fig.pdfC = 'none', info.txtC = 'none')
scores <- getScoreMN(pca)               # samples x components
getSummaryDF(pca)                       # R2X(cum) per component
# Tight pooled-QC clustering in the center = trustworthy run; QC scatter = analytical variance dominates
```

## Permutation-Validated PLS-DA / OPLS-DA

**Goal:** Decide whether group separation is real, not a geometry artifact, before reading any VIP or S-plot.

**Approach:** Fit with an explicit scaling, raise `permI` far above the default of 20, and read `pQ2`/`pR2Y` -- a model whose true Q2 sits inside the permutation cloud is indistinguishable from chance.

```r
library(ropls)
group <- factor(sample_info$group)
# OPLS-DA: 1 predictive + auto orthogonal; permI default 20 is too few for a reliable pQ2 -> >=1000
oplsda <- opls(t(feature_matrix), group, predI = 1, orthoI = NA,
               scaleC = 'pareto', permI = 1000, crossvalI = 7,
               fig.pdfC = 'none', info.txtC = 'none')
summ <- getSummaryDF(oplsda)            # R2X(cum), R2Y(cum), Q2(cum), pre, ort, pR2Y, pQ2
vip_pred <- getVipVn(oplsda)            # predictive VIP (Galindo-Prieto 2014); orthoL=TRUE for orthogonal
# Claim is licensed only if Q2 high AND pQ2 small. R2Y alone proves nothing.
```

PLS-DA is `orthoI = 0`. OPLS-DA has identical predictive power to PLS-DA -- it is a coordinate rotation, not a better model; the orthogonal block often encodes a confounder (inspect what correl