--- name: gaudi description: 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. category: analysis license: CC-BY-4.0 author: Heureka Labs version: 1.3.0 tags: [multi-omics, umap, hdbscan, clustering, shap] allowed-tools: Read, Write, Edit, Bash datasets: [] verified: date: 2026-08-29 against: Python route — umap-learn 0.5.12 / hdbscan 0.8.44 / xgboost 3.2.0 / shap 0.51.0 / scikit-learn 1.9.0 / NumPy 2.4.6 / pandas 3.0.5 / Python 3.11.15. R route — R 4.4.1 / gaudi 0.1.18 installed from GitHub / xgboost 3.2.1.1. Parameter defaults additionally cross-read against the gaudi 0.1.18 sources. executed: 6 unverified: 0 --- # 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 against gaudi 0.1.18 — the Python route by execution, the R route by executing the install, load and availability blocks and by reading the 0.1.18 sources for the parameter defaults.** 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`. ```r # install.packages("devtools") devtools::install_github("hirscheylab/gaudi") ``` Some functions need packages **attached**, not merely installed (see *Sharp edges*). The safe preamble: ```r library(gaudi) library(dplyr) # required by the survival helpers library(patchwork) # required by plot_gaudi_grid() ``` Check what is present before writing a script: ```bash 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**: ```r 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: ```r 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 ```r 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 six 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). ```r 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. ```r 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. ### `compute_features = TRUE` — the default — fails on xgboost R 3.x The default path dies before it returns, on any current R toolchain: ``` Error in `colnames<-`(`*tmp*`, value = c(colnames(X_train), "(Intercept)")) : attempt to set 'colnames' on an object with less than two dimensions Calls: gaudi -> xgboost_model -> -> colnames<- ``` Nothing in gaudi changed. `gaudi::xgboost_model` calls the pre-3.0 R interface (`xgboost(data =, label =, params =, verbose =)`); xgboost 3.x renamed `data` to `x` and `label` to `y` and dropped `params`, and its deprecation shims let training complete with the wrong arguments rather than stopping. `predict(mod, x, predcontrib = TRUE)` then hands back a plain numeric vector instead of an n x (p+1) matrix, and `SHAPforxgboost::shap.values` fails setting column names on a one-dimensional object. Everything downstream of the metagenes is unreachable while this holds: the `metagenes` slot, `plot_metagenes()`, `plot_gaudi_grid()` and `gaudi_enrichment()`. **`compute_features = FALSE` runs clean**, and gives you clusters and factors without the SHAP attribution. This is the failure a source read cannot find. The gaudi sources are unchanged and self-consistent; what moved was a dependency underneath them. Check it against your own installed xgboost rather than trusting either the package or this page. ### `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: ```r 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: ```r 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: ```r 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: ```r 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: ```r 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. ```r 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: ```r 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: ```r 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 ```r 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. ## Try it A self-contained check that this skill still works. No account, no key, no network beyond installing the packages. **Data** — generated inline, which is why `datasets:` is empty. The method needs a *known* partition to be checked against, and no public multi-omics dataset supplies one: three groups of twenty samples across two layers, with the signal in the first 30 features of each, is a case where the right answer is three clusters and anything else is a failure. The generator is pinned here so the figures below can be reproduced exactly rather than approximately. **What this does and does not exercise.** It runs the **Python** route from `references/python-implementation.md`, inlined so the block runs cold in an empty directory. The R package is not exercised — installing it needs an R toolchain — so the *Sharp edges* above remain source-confirmed rather than executed. **Run** — `pip install numpy pandas scikit-learn umap-learn hdbscan xgboost shap`, then: ```python import inspect import numpy as np, pandas as pd, umap, hdbscan from sklearn.metrics import silhouette_score, adjusted_rand_score def gaudi(omics, min_pts=None, n_components_layer=4, n_components_conc=2, n_neighbors=15, min_dist=0.01, compute_features=False, seed=42): """The implementation from references/python-implementation.md, metagenes off.""" common = sorted(set.intersection(*[set(df.index) for df in omics.values()])) if len(common) < 10: raise ValueError(f"only {len(common)} samples shared across layers") omics = {k: v.loc[common] for k, v in omics.items()} for k, v in omics.items(): if v.isna().any().any(): raise ValueError(f"NA in layer '{k}' — impute upstream") omics[k] = v.loc[:, v.std(axis=0, ddof=1) > 0] emb = [umap.UMAP(n_neighbors=n_neighbors, n_components=n_components_layer, min_dist=min_dist, metric="euclidean", random_state=seed).fit_transform(v.values) for v in omics.values()] integ = umap.UMAP(n_neighbors=n_neighbors, n_components=n_components_conc, min_dist=min_dist, metric="euclidean", random_state=seed).fit_transform(np.hstack(emb)) n = len(common) if min_pts is None: min_pts = int(np.floor(0.03 * n)) min_pts = max(min_pts, 2) lab = hdbscan.HDBSCAN(min_cluster_size=min_pts, min_samples=min_pts).fit_predict(integ) + 1 sil = float(silhouette_score(integ, lab)) if len(set(lab)) > 1 else 0.0 return dict(clusters=lab, silhouette_score=sil, n=n) def layers(shift, n_per=20, seed=0): """Three groups; signal in the first 30 features of each layer.""" rng = np.random.default_rng(seed) truth = np.repeat([0, 1, 2], n_per) ids = [f"S{i:02d}" for i in range(n_per * 3)] out = {} for name, p in (("expression", 300), ("methylation", 150)): X = rng.normal(size=(n_per * 3, p)) for g in range(3): X[truth == g, :30] += shift * (g - 1) out[name] = pd.DataFrame(X, index=ids, columns=[f"{name[:4]}{j}" for j in range(p)]) return out, truth k = lambda lab: len(set(lab) - {0}) # --- INVARIANTS: a failure here means the skill is wrong ---------------------- clean, truth = layers(shift=3.0) r = gaudi(clean, min_pts=5) assert k(r["clusters"]) == 3, k(r["clusters"]) # separable groups are recovered assert (r["clusters"] == 0).sum() == 0 # nothing left as noise assert adjusted_rand_score(truth, r["clusters"]) == 1.0 # exactly the true partition assert r["clusters"].min() >= 0 # R's convention: 0 = noise, never -1 print(f"clean (shift 3.0), min_pts=5 : k={k(r['clusters'])} noise=0 ARI=1.000 " f"sil={r['silhouette_score']:.3f}") # umap-learn's min_dist default is 10x uwot's, which GAUDI inherits. The # implementation passes 0.01 explicitly; if this ever fails, the porting note changed. assert inspect.signature(umap.UMAP).parameters["min_dist"].default == 0.1 print("umap-learn min_dist default : 0.1 (uwot's is 0.01 — pass it explicitly)") # Error paths: the Python route rejects what the R package silently drops. na = {n: d.copy() for n, d in clean.items()} na["expression"].iloc[3, 7] = np.nan try: gaudi(na, min_pts=5); raise AssertionError("NA should have raised") except ValueError as e: print("NA in a layer :", e) try: gaudi({n: d.iloc[:5] for n, d in clean.items()}, min_pts=2) raise AssertionError("5 samples should have raised") except ValueError as e: print("too few shared samples :", e) # --- OBSERVED, 2026-08-28: drift to investigate, not a bug ------------------- print("\nmin_pts x separation (truth k=3, n=60, default min_pts = 2):") print(f"{'shift':>6} {'default':>18} {'min_pts=5':>18} {'min_pts=10':>18}") for shift in (0.8, 3.0): om, tr = layers(shift=shift) cells = [] for mp in (None, 5, 10): rr = gaudi(om, min_pts=mp) cells.append(f"k={k(rr['clusters']):2d} ARI={adjusted_rand_score(tr, rr['clusters']):.2f}") print(f"{shift:>6} {cells[0]:>18} {cells[1]:>18} {cells[2]:>18}") # min_pts at or above the group size returns NOTHING, silently: no exception, # every sample labelled noise, and silhouette reported as 0.0. big = gaudi(clean, min_pts=20) print(f"\nmin_pts=20 (group size 20) : k={k(big['clusters'])} " f"noise={(big['clusters'] == 0).sum()} sil={big['silhouette_score']:.3f}") if k(big["clusters"]) == 0: assert big["silhouette_score"] == 0.0 # an empty result, reported as a number ``` **Expect** — the assertions above are invariants: on separable data the method returns the true partition, noise is labelled `0` and never `-1`, and both error paths raise rather than returning a quietly smaller result. The numbers below are **observed** on the versions in `verified:` and dated — a mismatch is drift to investigate, not a bug. UMAP prints `UserWarning: n_jobs value 1 overridden to 1 by setting random_state` once per embedding; that is expected, and it is the seed doing its job. ``` clean (shift 3.0), min_pts=5 : k=3 noise=0 ARI=1.000 sil=0.941 umap-learn min_dist default : 0.1 (uwot's is 0.01 — pass it explicitly) NA in a layer : NA in layer 'expression' — impute upstream too few shared samples : only 5 samples shared across layers min_pts x separation (truth k=3, n=60, default min_pts = 2): shift default min_pts=5 min_pts=10 0.8 k= 5 ARI=0.70 k= 3 ARI=0.75 k= 3 ARI=0.78 3.0 k= 3 ARI=1.00 k= 3 ARI=1.00 k= 3 ARI=1.00 min_pts=20 (group size 20) : k=0 noise=60 sil=0.000 ``` Two things in that table are the point of running it. The default `min_pts` is **not** safe merely because a clean dataset survives it — at `shift 3.0` it returns the truth, and at `shift 0.8` the same default splits three groups into five. Separation decides, and you do not know your separation in advance, which is why the sweep is not optional. And `min_pts=20`, at the group size, returns **every sample as noise with a silhouette of 0.0 and no error at all** — an empty result that reads like a computed one. **This collapse is specific to the Python route.** The same sweep run in R against gaudi 0.1.18 returns `k=3, noise=0` at every value from 5 to 20 on the same shape of data, so the exposure is shared but the size of the effect is not. Bound the sweep from above as well as below on the Python route, and read the noise count before the silhouette on either. ## 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 - Castellano-Escuder P., Zachman D.K., Han K., Hirschey M.D. GAUDI: interpretable multi-omics integration with UMAP embeddings and density-based clustering. *Nat Commun* 16, 5771 (2025). https://doi.org/10.1038/s41467-025-60822-1 - Package: https://github.com/hirscheylab/gaudi (GPL-3) - Benchmarks: https://github.com/hirscheylab/umap_multiomics_integration (MIT) - Archived code: https://doi.org/10.5281/zenodo.15442172 - McInnes L., Healy J., Melville J. UMAP. arXiv:1802.03426 - Campello R.J.G.B. et al. HDBSCAN. *PAKDD* (2013) - Lundberg S.M., Lee S.-I. A unified approach to interpreting model predictions. *NeurIPS* (2017)