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.
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
remotes::install_github("herbermann/mdclust")
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 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"
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"
)
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