← All skills

gaudi

analysis

Unsupervised multi-omics integration with GAUDI — two-stage UMAP, HDBSCAN clustering, and SHAP metagenes, in R or Python. Use when two or more omics layers share samples and you need sample clusters plus the features driving them.

GAUDI: multi-omics integration via UMAP + HDBSCAN

GAUDI (Group Aggregation via UMAP Data Integration) is an unsupervised, non-linear multi-omics integration method. It embeds each omics layer with UMAP independently, concatenates those embeddings, runs a second UMAP over the concatenation, clusters the result with HDBSCAN, and attributes features to the latent dimensions with SHAP values from XGBoost.

Reference implementation: the gaudi R package (GPL-3), Castellano-Escuder et al., Nat Commun 16:5771 (2025), doi:10.1038/s41467-025-60822-1.

Use it when you have two or more omics layers measured on the same samples and want sample clusters plus an interpretable account of which features drive them. It is a clustering and stratification method, not a supervised predictor and not an imputation method.

Everything documented below was verified by running gaudi 0.1.18 under R 4.3.2. The failure modes in Sharp edges are real and reproducible — read that section before running anything on real data. Several of them fail silently, and one of them silently deletes samples.

No R available, or Python-only? GAUDI is a composition of UMAP, HDBSCAN and SHAP — the R package is one implementation of the method, not the method itself. references/python-implementation.md states the algorithm language-agnostically and gives a tested Python version that reproduces the R package's clusters exactly on separable data. Do not install R just to run this. Everything in this document other than the R call syntax — the input contract, the parameter choices, the sharp edges, the reporting standards — applies to either route.

Why two stages of UMAP

The first UMAP normalises each layer's contribution: a 25,000-feature transcriptome and a 700-feature metabolome both become an n×4 embedding, so the integration is not dominated by whichever layer happens to be highest-dimensional. The second UMAP integrates the standardised representations. The paper validated this against nine alternative integration strategies (intersection, union, subtraction, joint matrix factorisation, UMAP on raw concatenated data) and selected concatenate-then-UMAP on cluster purity and silhouette.

HDBSCAN is used rather than k-means because the second UMAP's latent space is non-linear: HDBSCAN handles clusters of varying density and irregular shape, and does not require the cluster count in advance. The cost is a dedicated noise label (cluster 0) that you must handle deliberately.

Install and load

gaudi is not on CRAN or Bioconductor. If R is unavailable, skip this section entirely and use references/python-implementation.md.

# install.packages("devtools")
devtools::install_github("hirscheylab/gaudi")

Some functions need packages attached, not merely installed (see Sharp edges). The safe preamble:

library(gaudi)
library(dplyr)       # required by the survival helpers
library(patchwork)   # required by plot_gaudi_grid()

Check what is present before writing a script:

Rscript -e 'for (p in c("gaudi","uwot","dbscan","xgboost","SHAPforxgboost","dplyr","patchwork","sva","vegan","survminer","clusterProfiler","org.Hs.eg.db")) cat(sprintf("%-16s %s\n", p, ifelse(requireNamespace(p, quietly=TRUE), as.character(packageVersion(p)), "MISSING")))'

Input contract

A named list of matrices or data frames, samples in rows, features in columns, with row names set to sample IDs:

omics_list <- list(
  expression  = expr_mat,    # n samples x p1 genes
  methylation = meth_mat,    # n samples x p2 probes
  protein     = prot_mat     # n samples x p3 proteins
)

Name the list. Element names carry through to individual_factors and make on_omics indices legible — the metagenes slot is unnamed, so the list order is your only record of which layer is which.

Layers need not share features or dimensionality, but they must overlap in samples. GAUDI silently intersects on sample names and returns only the intersection, so a naming mismatch shows up as a quietly smaller result rather than an error. Check it yourself:

Reduce(intersect, lapply(omics_list, rownames))   # confirm before running

Pass samples_in_rows = FALSE if your matrices are features×samples. Row names are mandatory in either orientation — unnamed input fails with a confusing n_neighbors must be smaller than the dataset size.

Preprocess before GAUDI

GAUDI does no normalisation. Do it per layer, upstream:

  • Impute missing values. GAUDI cannot tolerate NA — see Sharp edges.
  • Scale/standardise per layer (z-score is the usual choice).
  • Apply the layer-appropriate transform first: log-CPM or VST for counts, M-values for methylation betas, log for intensity-based metabolomics.
  • Drop features that are constant or near-constant.

Running it

res <- gaudi(
  omics_list,
  umap_params      = list(n_neighbors = 15, n_components = 4),  # per-layer
  umap_params_conc = list(n_neighbors = 15, n_components = 2),  # integrated
  min_pts          = NULL,       # NULL -> max(floor(0.03 * n), 2)  <- SET THIS
  xgboost_params   = list(lambda = 0, eta = 0.5, gamma = 50,
                          max_depth = 10, subsample = 0.95),
  compute_features = TRUE,
  combine_omics    = FALSE,
  samples_in_rows  = TRUE,
  reassign_cluster_zero = FALSE,
  method           = "xgboost"   # or "rf"
)

metric = "euclidean" and min_dist = 0.01 are not arguments — they are uwot::umap defaults that GAUDI inherits. To change them, add them to umap_params.

Published defaults and their stated rationale:

parameter default why
n_neighbors 15 uwot default; balances local/global structure
n_components (stage 1) 4 enough for primary variance components per layer
n_components (stage 2) 2 visualisation and interpretation
min_dist 0.01 tight local packing, favours compact clusters
lambda 0 no regularisation — sparse biological signal
eta 0.5 moderate learning rate
gamma 50 sparse trees; few features retained
max_depth 10 admits feature interactions
subsample 0.95 stable feature selection across runs

Sizing guidance from the package: reduce n_neighbors below 100 samples, consider raising it above 1000. Scaling is roughly linear; 1000 samples × 25,000 features × 5 layers runs in under 10 minutes on standard hardware.

Result object

GAUDIObject, an S4 object with five slots:

slot contents
@factors data.frame, n × 3: UMAP1, UMAP2, clust; row names = sample IDs
@clusters numeric vector of cluster labels; 0 = HDBSCAN noise
@silhouette_score mean silhouette over the 2-D embedding; 0 if one cluster
@individual_factors named list of per-layer n×4 embeddings
@metagenes unnamed list, one data.frame per layer: contrib1, contrib2

metagenes[[i]] holds mean absolute SHAP contributions, so values are non-negative — magnitude only, no direction. Rows are sorted by |contrib1| descending. Because gamma = 50 forces sparse trees, most features are exactly 0; that is the intended behaviour, not a failure. To recover direction, compare the feature across clusters directly (e.g. limma).

res@silhouette_score
table(res@clusters)                       # inspect the size of cluster 0
head(res@metagenes[[1]], 20)              # top features for layer 1
nonzero <- sum(res@metagenes[[1]]$contrib1 > 0)

Consensus variant

c_gaudi() runs gaudi() n_max times, picks the highest-silhouette run as a reference, Procrustes-aligns the others onto it (vegan::procrustes), takes the per-sample median embedding, and clusters that. Use it when UMAP run-to-run variability matters — reporting a single stochastic embedding as the result is not defensible.

res_c <- c_gaudi(omics_list, n_max = 10, min_pts = 12)

Two differences to know: c_gaudi defaults to different XGBoost parameters (eta = 1, gamma = 100 — sparser than gaudi()), and it returns individual_factors empty. Pass xgboost_params explicitly if you want metagenes comparable to gaudi().

Sharp edges

These are verified defects and traps in gaudi 0.1.18, not hypotheticals.

NA silently deletes samples, then crashes

align_omics() computes apply(x, 2, sd) > 0 on the features×samples matrix and keeps columns passing the test. A single NA anywhere in a sample's column makes its sd return NA, the test yields NA, and that sample is dropped. Because layers are filtered independently, they end up with different sample counts and the run dies in dplyr::bind_cols with:

Can't recycle `..1` (size 59) to match `..2` (size 60).

That message names neither NA nor the sample. Impute before calling GAUDI and assert completeness:

stopifnot(!any(sapply(omics_list, function(x) anyNA(x))))

For >20% missingness, drop the affected samples or features rather than imputing.

The zero-SD filter removes constant samples, not constant features

The same line intends to drop zero-variance features but, given the matrix orientation at that point, tests samples. Consequences: a constant/flatlined sample is dropped without warning, and genuinely constant features survive into the model. So filter zero-variance features yourself upstream, and verify the returned sample count:

stopifnot(nrow(res@factors) == length(Reduce(intersect, lapply(omics_list, rownames))))

The default min_pts fragments small datasets

min_pts defaults to max(floor(0.03 * n), 2). Below ~100 samples this collapses to 2 or 3, and HDBSCAN returns a shower of doubleton clusters. On a 60-sample, 3-group synthetic set with clean separation:

min_pts clusters found noise silhouette
2 (the default at n=60) 19 4 0.33
5 3 0 0.96
10 3 0 0.96
18 3 0 0.94

The truth was 3. Always set min_pts explicitly below ~150 samples, and sweep it rather than accepting the first answer:

for (mp in c(5, 8, 10, 15, 20)) {
  r <- gaudi(omics_list, min_pts = mp, compute_features = FALSE)
  cat(sprintf("min_pts=%2d  k=%2d  noise=%3d  sil=%.3f\n",
      mp, length(setdiff(unique(r@clusters), 0)), sum(r@clusters == 0),
      r@silhouette_score))
}

Report the sweep, not just the chosen value. A cluster count that swings wildly across min_pts is a finding about the data, and silhouette alone will not tell you the count is wrong — it is computed on the same embedding that produced the clusters, so it rewards the geometry UMAP already imposed. Anchor cluster counts to something external (survival, known subtypes, held-out annotation) before claiming they are biological.

gaudi_enrichment() only works on factor 1

Metagene rows are sorted by |contrib1|, and clusterProfiler::gseGO requires a decreasing-sorted input. So on_factor = 2 fails with geneList should be a decreasing sorted vector. Sort it yourself for dimension 2:

gl <- res@metagenes[[1]][["contrib2"]]
names(gl) <- rownames(res@metagenes[[1]])
gl <- sort(gl[gl != 0], decreasing = TRUE)
gse <- clusterProfiler::gseGO(gl, ont = "ALL", OrgDb = org.Hs.eg.db,
                              keyType = "SYMBOL", minGSSize = 2, maxGSSize = 200)

Also: gaudi_enrichment() hardcodes pvalueCutoff = 1 and pAdjustMethod = "none", so it returns everything unfiltered and uncorrected. You must correct and filter yourself:

e <- gaudi_enrichment(res, on_omics = 1, organism = org.Hs.eg.db)
e$padj <- p.adjust(e$pvalue, method = "BH")
e <- e[e$padj < 0.05, ]

Never report rows straight out of that function as significant. Note also that SHAP scores are all non-negative, so GSEA is running on a one-sided statistic — fgsea will warn about scoreType. Treat the result as ranked-magnitude enrichment, not directional enrichment.

Row names of the chosen layer must be gene symbols for this to work.

Survival helpers need data frames and dplyr attached

get_pairwise_survival_data() and plot_survival() are documented as accepting named vectors. They do not — the vector path hits an unquoted := and dies with could not find function ":=". Pass data frames with a time / censor column and sample row names. Both functions also call bare mutate() / across(), so library(dplyr) must be attached.

library(dplyr)
time_df   <- data.frame(time   = clin$os_days,   row.names = clin$sample_id)
censor_df <- data.frame(censor = clin$os_event,  row.names = clin$sample_id)

pairs <- get_pairwise_survival_data(res, time_df, censor_df)
head(pairs)   # clusters, pval, distance, quality

plot_survival(res, time_df, censor_df, pairs[["clusters"]][[1]])

quality is 0.8 * normalised(-log p) + 0.2 * scaled centroid distance — an uninterpretable triage heuristic for ranking which cluster pair to examine, not a statistic. The pval column is a raw pairwise log-rank p-value across all choose(k, 2) comparisons; correct it before reporting:

pairs$padj <- p.adjust(pairs$pval, method = "BH")

plot_survival() drops cluster 0 and any NA rows automatically.

plot_gaudi_grid() needs patchwork attached

It composes panels with | and /. Without library(patchwork) it fails with non-numeric argument to binary operator. plot_factors() and plot_metagenes() are unaffected.

combine_omics = TRUE shifts every metagene index

With combine_omics = TRUE, GAUDI ComBat-corrects the concatenated layers (treating layer identity as batch) and prepends that combined result to metagenes. For a two-layer input, length(metagenes) becomes 3: [[1]] = combined, [[2]] = layer 1, [[3]] = layer 2. Every on_omics index shifts by one. Verify by row count:

sapply(res@metagenes, nrow)

Cluster 0 is noise, and dropping it is a decision

drop_clusters(res, clusters = 0) removes noise samples; reassign_cluster_zero = TRUE instead assigns each to its nearest neighbour's cluster. Both change downstream n. State which you used and why — silently absorbing noise into clusters inflates apparent separation. Report the noise fraction either way.

Note that reassign_cluster_zero = TRUE reorders @factors rows (noise samples move to the end), so re-match by row name rather than position when joining metadata.

A defensible run

library(gaudi); library(dplyr); library(patchwork)
set.seed(42)   # UMAP and XGBoost are both stochastic

stopifnot(!any(sapply(omics_list, anyNA)))
common <- Reduce(intersect, lapply(omics_list, rownames))
message(length(common), " samples shared across ", length(omics_list), " layers")

res <- gaudi(omics_list, min_pts = 12, samples_in_rows = TRUE)

cat("k =", length(setdiff(unique(res@clusters), 0)),
    "| noise =", sum(res@clusters == 0),
    "| silhouette =", round(res@silhouette_score, 3), "\n")
stopifnot(nrow(res@factors) == length(common))

ggplot2::ggsave("results/gaudi_embedding.png", plot_factors(res, label_size = 4),
                width = 7, height = 6, dpi = 300)
ggplot2::ggsave("results/gaudi_grid.png", plot_gaudi_grid(res, top_features = 15),
                width = 14, height = 7, dpi = 300)

clusters <- tibble::rownames_to_column(res@factors, "sample_id")
readr::write_csv(clusters, "results/gaudi_clusters.csv")

Record for provenance — min_pts (and the sweep), the seed, the per-layer preprocessing, gaudi version (res@gaudiVersion), noise fraction, silhouette, and how cluster 0 was handled.

For a guarded end-to-end version of this — the same run with every failure mode in Sharp edges turned into a loud error, plus a min_pts sweep and a provenance file — see references/analysis-template.md. Copy the block into gaudi_analysis.R, fill in the CONFIG section, and run it with Rscript.

Reporting

State plainly:

  • number of layers, features per layer, and the final sample n after intersection
  • preprocessing and imputation per layer
  • min_pts and how it was chosen; the cluster-count sweep
  • clusters found, noise fraction, treatment of cluster 0
  • silhouette — and that it is computed on the same embedding it evaluates
  • random seed; whether gaudi() or c_gaudi() was used
  • for metagenes: that SHAP values are unsigned importance
  • multiple-testing correction for pairwise survival and enrichment

Do not claim clusters are biological subtypes on silhouette alone. Do not read UMAP inter-cluster distances as quantitative — UMAP preserves local neighbourhoods, not global geometry, so the layout is topology, not a metric space. Cluster assignments are the output; the coordinates are a picture of them.

Limits

  • Unsupervised clustering, not classification or prediction. No predict(); new samples require re-running the whole integration.
  • Stochastic. Set a seed; prefer c_gaudi() when stability matters.
  • Complete data required. No NA tolerance — impute upstream.
  • Requires matched samples across layers. For unmatched single-cell layers, the paper matches cells first (PCA + greedy matching, validated with Procrustes) before GAUDI; that matching is a separate analysis with its own assumptions.
  • Unsigned feature attribution. Magnitude only.
  • Silhouette is internal. It cannot validate that the cluster count is biologically real.

Reference files

  • references/api-reference.md — every exported signature, the full default set for gaudi() and c_gaudi(), the GAUDIObject slot table, and a behaviour matrix of what was observed empirically. Read it when you need an argument you cannot find above, or to confirm a slot's shape before indexing.
  • references/analysis-template.md — a complete guarded analysis script. Copy the fenced block to gaudi_analysis.R and edit the CONFIG section.
  • references/python-implementation.md — the method stated language-agnostically, plus a tested Python implementation for when R is unavailable. Includes the porting traps that make a naive translation silently wrong (min_dist defaults differ 10×; the noise label differs) and a verified account of what reproduces across implementations and what does not.

Sources