mdclust

mdclust

mdclust is a standalone implementation of Maximum Diversity Clustering (Algorithm 2 from Harmony), providing a reusable clustering engine with rich output and complete control over the algorithm.

Features

• Standalone implementation of Maximum Diversity Clustering (Algorithm 2 from Harmony)

• Soft spherical k-means with a batch-diversity penalty

• k-means++ initialization

• Rich optimization diagnostics (O, E, logOE, objective decomposition, …)

• Algorithm parameters fully configurable for more control and use in complex pipelines

• Native C++ implementation via RcppArmadillo

Installation

remotes::install_github("herbermann/mdclust")

Minimal example

Simulated data

We use simulated data to showcase the functionality.

# --- Simulate signal // groups A & B ---
n <- 600

signal <- rep(c("A", "B"), each = n / 2)
Z <- matrix(0, nrow = n, ncol = 2)
Z[signal == "A", ] <- cbind(rnorm(n / 2, -0.8, 0.8), rnorm(n / 2, 0.0, 0.5))
Z[signal == "B", ] <- cbind(rnorm(n / 2, 0.8, 0.8), rnorm(n / 2, 0.0, 0.5))

# --- Assign batches independently ---
batch <- sample(c("Batch 1", "Batch 2"), n, replace = TRUE)

# --- Add strong batch effect ---
shift <- c(0., 2.)
Z[batch == "Batch 2", ] <- sweep(Z[batch == "Batch 2", , drop = FALSE], 2, shift, "+")

Running mdclust

Running mdclust is as simple as

library(mdclust)

fit <- mdclust(
  t(Z),
  Y_init = matrix(rnorm(2 * 2), 2, 2),
  batch = batch,
  K = 2
)

fit

<mdclust>

 Elements        : 600
 Clusters        : 2
 Embedding dims  : 2
 Converged       : yes
 Iterations      : 11
 Final objective : 251.015

Components:
  R           soft assignments
  cluster     hard assignments
  Y           centroids
  objective   optimization trace
  O, E, logOE batch diagnostics

Here, we provided an explicit initial centroid position Y_init. Other options are "kmeans", and "kmeans++" respectively. If none is provided, mdclust makes an optimal initial guess using (spherical) kmeans++. More parameters of the algorithm can be chosen freely, for more details check also the vignettes.

The returned object provides rich access to internals, and the optimization details.

names(fit)

[1] "R"            "Y"            "O"            "E"            "logOE"       
[6] "objective"    "converged"    "cluster"      "batch_levels"

Relationship to Harmony

Maximum Diversity Clustering is one component of the Harmony integration algorithm. Harmony exposes this clustering only as an internal step of the integration procedure. While the resulting clustering state (R, Y, O, E, …) can be inspected by requesting the full Harmony object, the clustering itself is not available as a standalone, fully configurable algorithm. In particular, Harmony does not expose initialization of the cluster centroids or alternative optimization strategies through its public API. Full control over algorithm parameters as well as the full output might be desirable for some standalone uses or integration in larger pipelines. Offering this is the primary intent of this package.

A roughly equivalent call (up to potential discrepancies in the way Y_init is chosen) in harmony would be

set.seed(1)
library(harmony)

opts <- harmony::harmony_options(
  max.iter.cluster = 200,
  epsilon.cluster = 1e-5
)

harmony_object <- harmony::RunHarmony(
  data_mat = t(Z),
  meta_data = data.frame(batch = batch),
  vars_use = "batch",
  nclust = 2,
  theta = 2,
  sigma = 0.1,
  max_iter = 0,
  return_object = TRUE,
  .options = opts
)

Initializing centroids

A direct comparison shows mdclust works as intended.

fit <- mdclust(
  t(Z),
  batch = batch,
  K = 2
)


par(mfrow = c(1, 3), mar = c(4, 4, 3, 1))

plot(
  Z,
  col = factor(signal),
  pch = 16,
  asp = 1,
  main = "Ground truth",
  xlab = "Dim 1",
  ylab = "Dim 2"
)

harmony_cluster <- max.col(t(harmony_object$getR()))

plot(
  Z,
  col = factor(harmony_cluster),
  pch = 16,
  asp = 1,
  main = "Harmony",
  xlab = "Dim 1",
  ylab = "Dim 2"
)
plot(
  Z,
  col = factor(fit$cluster),
  pch = 16,
  asp = 1,
  main = "mdclust",
  xlab = "Dim 1",
  ylab = "Dim 2"
)

Citation

Korsunsky et al. (2019). Fast, sensitive and accurate integration of single-cell data with Harmony. Nature Methods, 16, 1289–1296. https://doi.org/10.1038/s41592-019-0619-0