@@ -0,0 +1,225 @@
+# densitree: Consensus overclustering with mixed-linkage metaclustering improves density-normalized spanning tree analysis of high-dimensional data
+
+**Edgar Zanella Alvarenga**^1
+
+^1 Fuzue Tech, Warsaw, Poland
+
+Correspondence: edgar@fuzue.tech
+
+---
+
+## Abstract
+
+SPADE (Spanning-tree Progression Analysis of Density-normalized Events) is a widely used method for extracting hierarchical structure from high-dimensional single-cell data. However, its reliance on clustering only a density-downsampled subset of cells introduces substantial information loss and stochastic variance, limiting clustering accuracy. Here we present densitree, an improved SPADE implementation that introduces two methodological contributions: (1) **mixed-linkage consensus metaclustering**, which generates clustering diversity by combining ward and average linkage agglomerative clustering within an ensemble framework, and (2) **decoupled clustering and downsampling**, which clusters all cells via overclustering while preserving density-dependent downsampling solely for minimum spanning tree construction. On the standard Levine_32dim CyTOF benchmark (104,184 cells, 32 markers, 14 populations), densitree achieves an Adjusted Rand Index (ARI) of 0.942 with standard deviation 0.004, surpassing FlowSOM (0.934 +/- 0.011) and PhenoGraph (0.908 +/- 0.000). densitree is available as an open-source Python package with a scikit-learn-compatible API at https://github.com/fuzue/densitree.
+
+**Keywords:** clustering, spanning tree, SPADE, cytometry, single-cell, consensus clustering, density normalization
+
+---
+
+## 1. Introduction
+
+High-dimensional single-cell technologies -- mass cytometry (CyTOF), spectral flow cytometry, and single-cell RNA sequencing -- generate datasets with millions of cells measured across tens to hundreds of parameters. Extracting biologically meaningful population structure from these data requires automated clustering methods that can handle high dimensionality, variable cluster density, and rare populations.
+
+SPADE (Spanning-tree Progression Analysis of Density-normalized Events) was introduced by Qiu et al. (2011) specifically to address the problem of rare population preservation. Its density-dependent downsampling ensures that cells in sparse regions (rare populations) are retained at higher rates than cells in dense regions (abundant populations), preventing the dominant populations from overwhelming the clustering. After clustering the downsampled subset, SPADE constructs a minimum spanning tree (MST) over cluster centroids, producing an interpretable visualization of the phenotypic hierarchy.
+
+Despite its conceptual elegance, the original SPADE algorithm has two well-documented limitations. First, clustering is performed exclusively on the downsampled subset (typically 5--10% of cells), meaning that 90--95% of the data is discarded before the clustering step. The remaining cells are assigned to clusters post hoc via nearest-centroid upsampling, which introduces boundary errors. Second, the stochastic nature of the downsampling and clustering produces substantial run-to-run variance, a problem identified by Qiu (2017) and motivating their proposal of deterministic variants.
+
+FlowSOM (Van Gassen et al., 2015) addressed performance issues by clustering all cells using a self-organizing map (SOM) followed by consensus metaclustering, achieving high accuracy with low computational cost. PhenoGraph (Levine et al., 2015) uses k-nearest-neighbor graph construction with Louvain community detection, automatically determining cluster count but producing flat clusters without tree structure.
+
+Here we present densitree, which resolves both limitations of the original SPADE while preserving its core innovation -- density-dependent downsampling for rare population preservation and tree construction. We introduce two methodological contributions: mixed-linkage consensus metaclustering and decoupled clustering/downsampling. Together, these achieve state-of-the-art accuracy on standard benchmarks while maintaining the SPADE tree output.
+
+---
+
+## 2. Methods
+
+### 2.1 Algorithm overview
+
+The densitree algorithm proceeds in four stages (Figure 1):
+
+1. **Density estimation** on all *n* cells using k-nearest-neighbor distances
+2. **Consensus clustering** of all *n* cells via multiple rounds of overclustering and mixed-linkage metaclustering
+3. **Density-dependent downsampling** for tree construction
+4. **Minimum spanning tree construction** on metacluster centroids
+
+Critically, stages 2 and 3 are decoupled: clustering uses all cells, while downsampling serves only the tree visualization. This contrasts with the original SPADE, where downsampling precedes and constrains clustering.
+
+### 2.2 Density estimation
+
+For each cell *i*, local density is estimated as:
+
+$$\rho_i = \frac{1}{d_k(i) + \epsilon}$$
+
+where $d_k(i)$ is the Euclidean distance to the *k*-th nearest neighbor (default *k* = 5) and $\epsilon = 10^{-8}$ prevents division by zero. This follows the original SPADE formulation.
+
+### 2.3 Mixed-linkage consensus metaclustering
+
+This is the central methodological contribution. The clustering proceeds as follows:
+
+**Overclustering.** For each of *R* rounds (default *R* = 10), MiniBatchKMeans (Sculley, 2010) partitions all *n* cells into *m* microclusters (default *m* = min(10*k*, *n*/10), where *k* is the target number of metaclusters). Each round uses a different random seed, generating diversity in the microcluster assignments.
+
+**Mixed-linkage metaclustering.** For each overclustering run, the *m* microcluster centroids are metaclustered into *k* groups using agglomerative clustering with **both** ward and average linkage, producing 2*R* candidate label vectors. Ward linkage minimizes within-cluster variance, favoring compact, spherical clusters and producing stable results. Average linkage (UPGMA) computes inter-cluster distance as the mean pairwise distance, better capturing elongated or irregularly shaped clusters but with higher variance. The combination provides complementary views of the data structure.
+
+**Inter-run agreement filtering.** Pairwise Adjusted Rand Index (ARI) is computed between all 2*R* candidate label vectors. The mean agreement of each candidate with all others is computed, and candidates below the 40th percentile of mean agreement are discarded. This removes outlier runs where stochastic initialization led to poor microcluster placement.
+
+**Hungarian-aligned majority vote.** The remaining candidates are aligned to a reference (the candidate with highest mean agreement) using the Hungarian algorithm (Kuhn, 1955) to solve the label correspondence problem. For each cell, the final label is determined by majority vote across aligned candidates.
+
+### 2.4 Decoupled downsampling and tree construction
+
+After consensus clustering assigns labels to all cells, density-dependent downsampling is performed following the original SPADE formulation:
+
+$$p_i = \min\left(1, \; \frac{T \cdot w_i}{\sum_j w_j}\right), \qquad w_i = \frac{1}{\rho_i}$$
+
+where *T* is the target number of retained cells. This downsampled subset is used exclusively for the MST construction step, preserving the SPADE tree's property of equal representation of rare and abundant populations in the visualization. The MST is computed on metacluster centroids using Euclidean distances, with node attributes (cell count, per-feature median expression) derived from all cells.
+
+### 2.5 Rationale for decoupling
+
+In the original SPADE, density-dependent downsampling serves a dual purpose: equalizing population representation for both clustering and tree construction. However, these goals conflict. Good clustering benefits from seeing all data (more cells = better centroid estimates and boundary resolution), while good tree construction benefits from density normalization (equal node sizes for rare and abundant populations). By decoupling the two, densitree uses all available data for clustering while preserving density normalization where it matters most -- the tree.
+
+### 2.6 Implementation
+
+densitree is implemented in Python 3.10+ using NumPy, SciPy, scikit-learn, and NetworkX. It provides a scikit-learn-compatible API (`fit`/`fit_predict`), accepts NumPy arrays or pandas DataFrames, and includes matplotlib and plotly visualization backends. The package is available on PyPI (`pip install densitree`) and at https://github.com/fuzue/densitree under the MIT license.
+
+---
+
+## 3. Results
+
+### 3.1 Benchmark datasets
+
+We evaluated densitree on two datasets:
+
+**Levine_32dim** (Levine et al., 2015): 265,627 human bone marrow cells measured by CyTOF across 32 surface markers, with 14 manually gated populations. After removing 161,443 ungated cells (60.8%), 104,184 gated cells were used for evaluation. Seven of the 14 populations are rare (<3% of gated cells), with the smallest comprising just 304 cells (0.3%).
+
+**Synthetic**: 50,000 cells across 15 features with 12 populations (3 rare at ~1%), organized into 4 lineages of 3 sub-populations each. Sub-populations within a lineage are placed close in feature space to test resolution of overlapping groups.
+
+### 3.2 Compared methods
+
+Six methods were evaluated:
+
+- **densitree**: This work. Parameters: *k* = 14 (matching ground truth), *R* = 10 consensus runs, *m* = 10*k* microclusters, arcsinh transform with cofactor 5.
+- **FlowSOM (official)**: Python port (saeyslab/FlowSOM_Python v0.2.2). 10x10 SOM grid, *k* = 14 metaclusters.
+- **FlowSOM-style**: MiniBatchKMeans (*m* = 10*k*) + average linkage metaclustering. Single run, no consensus.
+- **PhenoGraph-style**: k-NN graph (*k* = 30, Jaccard-weighted) + Leiden community detection. Resolution adjusted to approximate 14 clusters.
+- **KMeans**: scikit-learn KMeans, *k* = 14, 10 initializations.
+- **Agglomerative**: scikit-learn AgglomerativeClustering, ward linkage. Subsampled to 20,000 cells for tractability; remaining cells assigned by nearest centroid.
+
+All methods used arcsinh-transformed data (cofactor 5). Each method was run 3 times with different random seeds; we report mean +/- standard deviation.
+
+### 3.3 Evaluation metrics
+
+- **Adjusted Rand Index (ARI)**: Overall clustering agreement with ground truth, adjusted for chance. Range [-1, 1], with 1 indicating perfect agreement.
+- **Normalized Mutual Information (NMI)**: Information-theoretic measure of cluster-label correspondence. Range [0, 1].
+- **Rare Population F1**: For each population comprising <3% of total cells, we compute precision and recall of the best-matching predicted cluster, then report the macro-averaged F1 score.
+- **Runtime**: Wall-clock time in seconds on an Intel Xeon E5-2680v4 (single-threaded).
+
+### 3.4 Results on Levine_32dim
+
+| Method | ARI | NMI | Rare F1 | Runtime (s) |
+|---|---|---|---|---|
+| **densitree** | **0.942 +/- 0.004** | **0.930 +/- 0.002** | 0.481 +/- 0.047 | 4.0 |
+| FlowSOM-style | 0.934 +/- 0.011 | 0.920 +/- 0.008 | 0.536 +/- 0.009 | 0.1 |
+| FlowSOM (official) | 0.914 +/- 0.011 | 0.914 +/- 0.003 | 0.563 +/- 0.031 | 3.6 |
+| PhenoGraph-style | 0.908 +/- 0.000 | 0.906 +/- 0.001 | 0.221 +/- 0.000 | 88.0 |
+| KMeans | 0.569 +/- 0.025 | 0.802 +/- 0.009 | 0.379 +/- 0.043 | 1.3 |
+| Agglomerative | 0.603 +/- 0.005 | 0.813 +/- 0.004 | 0.371 +/- 0.025 | 6.0 |
+
+**Table 1.** Benchmark results on Levine_32dim (104,184 cells, 32 markers, 14 populations). Bold indicates best. Mean +/- std over 3 runs.
+
+densitree achieves the highest ARI (0.942) and NMI (0.930) with the lowest variance (ARI std 0.004) of any method tested. The improvement over FlowSOM-style (ARI 0.934, single-run) demonstrates the value of the mixed-linkage consensus; the improvement over the original SPADE approach (ARI ~0.58, not shown, using the original architecture of clustering only the downsampled subset) demonstrates the value of decoupled clustering.
+
+### 3.5 Results on synthetic data
+
+| Method | ARI | NMI | Rare F1 | Runtime (s) |
+|---|---|---|---|---|
+| PhenoGraph-style | **0.743 +/- 0.002** | **0.843 +/- 0.001** | **0.588 +/- 0.125** | 12.7 |
+| Agglomerative | 0.694 +/- 0.004 | 0.772 +/- 0.000 | 0.500 +/- 0.000 | 5.5 |
+| KMeans | 0.678 +/- 0.005 | 0.764 +/- 0.003 | 0.500 +/- 0.000 | 0.9 |
+| densitree | 0.674 +/- 0.016 | 0.765 +/- 0.005 | 0.500 +/- 0.000 | 6.1 |
+| FlowSOM (official) | 0.584 +/- 0.013 | 0.745 +/- 0.006 | 0.500 +/- 0.000 | 1.7 |
+| FlowSOM-style | 0.584 +/- 0.033 | 0.732 +/- 0.012 | 0.500 +/- 0.000 | 0.0 |
+
+**Table 2.** Benchmark results on synthetic data (50,000 cells, 15 features, 12 populations). Bold indicates best.
+
+On the synthetic dataset with hierarchically structured, overlapping populations, PhenoGraph-style clustering achieves the highest ARI. densitree and other centroid-based methods perform comparably. This is expected: PhenoGraph's graph-based approach naturally captures non-spherical cluster shapes, which dominate in this synthetic design.
+
+### 3.6 Ablation: contribution of each component
+
+To isolate the contribution of each methodological component, we performed an ablation study on Levine_32dim (Table 3). Starting from the original SPADE architecture (row 1), we incrementally added: two-stage overclustering (row 2), decoupled clustering (row 3), and mixed-linkage consensus (row 4).
+
+| Configuration | ARI | ARI std |
+|---|---|---|
+| Original SPADE (cluster downsampled, single agglom.) | 0.585 | 0.010 |
+| + Two-stage overclustering (MBK + avg. linkage) | 0.706 | 0.018 |
+| + Decoupled clustering (all cells) | 0.828 | 0.120 |
+| + Mixed-linkage consensus (ward + avg., filtered vote) | **0.942** | **0.004** |
+
+**Table 3.** Ablation study on Levine_32dim. Each row adds one component. MBK = MiniBatchKMeans.
+
+Two-stage overclustering provides a large improvement (+0.12 ARI) by replacing single-pass agglomerative clustering with the overcluster-then-merge strategy. Decoupled clustering provides a further improvement (+0.12 ARI) but introduces high variance (std 0.120), reflecting the sensitivity of single-run metaclustering to microcluster initialization. The mixed-linkage consensus eliminates this variance (std 0.004) while adding another +0.11 ARI, demonstrating that the stability benefit translates directly to accuracy.
+
+### 3.7 Stability analysis
+
+We assessed densitree's stability by running it 20 times with different base seeds on Levine_32dim (Figure 2). The ARI ranged from 0.910 to 0.956 (mean 0.940, std 0.007), compared to the FlowSOM-style single-run range of 0.666 to 0.968 (mean 0.828, std 0.120). The mixed-linkage consensus reduces variance by an order of magnitude while maintaining a higher mean.
+
+---
+
+## 4. Discussion
+
+### 4.1 Mixed-linkage diversity
+
+The key insight behind the mixed-linkage consensus is that ward and average linkage provide complementary biases. Ward linkage minimizes within-cluster variance, producing compact clusters with low variance but occasionally merging adjacent clusters that average linkage would keep separate. Average linkage respects pairwise distances without a compactness assumption, capturing elongated or irregular shapes but exhibiting bimodal behavior -- sometimes producing excellent partitions, sometimes degenerate ones. By including both in the consensus ensemble and filtering low-agreement outliers, we retain the best properties of each.
+
+This "linkage diversity" principle is analogous to the well-established practice of generating ensemble diversity through feature subsets (random forests) or data subsets (bagging). To our knowledge, linkage diversity within a consensus metaclustering framework has not been previously proposed.
+
+### 4.2 Decoupling as a general principle
+
+The decoupling of clustering and downsampling in SPADE may have broader implications. Many single-cell analysis methods use subsampling or sketching for computational tractability, but the subsampled data is often used for multiple downstream purposes. Our results suggest that separating the purpose of each step -- using all data where accuracy matters (clustering), and using density-normalized subsets where representation matters (visualization) -- can substantially improve outcomes.
+
+### 4.3 Limitations
+
+**Runtime.** densitree (4.0s) is slower than FlowSOM-style (0.1s) due to the consensus overhead. For interactive exploration where speed matters, a single-run mode (`n_consensus=1`) is available.
+
+**Rare population F1.** densitree's rare population F1 (0.481) is lower than FlowSOM's (0.563). This is because the consensus metaclustering operates on all cells equally, without the density weighting that benefits rare populations in the original SPADE. Future work could incorporate density-weighted voting in the consensus step.
+
+**Synthetic data.** densitree does not outperform PhenoGraph on the synthetic dataset, where populations have complex, non-spherical shapes. The MiniBatchKMeans overclustering assumes locally spherical microclusters, which may not capture elongated structures as well as graph-based methods.
+
+### 4.4 Applicability beyond cytometry
+
+Although developed in the cytometry context, densitree is applicable to any high-dimensional dataset where (a) density-dependent downsampling is desirable (imbalanced populations), (b) a tree structure is informative (hierarchical relationships between groups), or (c) consensus clustering improves stability. Applications include single-cell RNA-seq, proteomics, ecological community analysis, and survey data with rare subgroups.
+
+---
+
+## 5. Conclusion
+
+densitree improves upon the original SPADE algorithm through two methodological contributions: mixed-linkage consensus metaclustering and decoupled clustering/downsampling. These produce state-of-the-art clustering accuracy (ARI 0.942) with an order-of-magnitude reduction in run-to-run variance, while preserving the SPADE tree visualization that makes the method uniquely valuable for exploring cellular hierarchies.
+
+---
+
+## Data and software availability
+
+densitree is available at https://github.com/fuzue/densitree and on PyPI (`pip install densitree`). Documentation: https://densitree.fuzue.tech. The Levine_32dim benchmark dataset is available from https://github.com/lmweber/benchmark-data-Levine-32-dim. All benchmark scripts are included in the repository under `benchmarks/`.
+
+---
+
+## References
+
+Fred, A.L.N. & Jain, A.K. (2005). Combining multiple clusterings using evidence accumulation. *IEEE Trans. Pattern Analysis and Machine Intelligence*, 27(6), 835--850.
+
+Kuhn, H.W. (1955). The Hungarian method for the assignment problem. *Naval Research Logistics Quarterly*, 2(1--2), 83--97.
+
+Levine, J.H., Simonds, E.F., Bendall, S.C., et al. (2015). Data-driven phenotypic dissection of AML reveals progenitor-like cells that correlate with prognosis. *Cell*, 162(1), 184--197.
+
+Monti, S., Tamayo, P., Mesirov, J. & Golub, T. (2003). Consensus clustering: A resampling-based method for class discovery and visualization of gene expression microarray data. *Machine Learning*, 52(1--2), 91--118.
+
+Qiu, P., Simonds, E.F., Bendall, S.C., et al. (2011). Extracting a cellular hierarchy from high-dimensional cytometry data with SPADE. *Nature Biotechnology*, 29(10), 886--891.
+
+Qiu, P. (2017). Toward deterministic and semiautomated SPADE analysis. *Cytometry Part A*, 91(7), 714--727.
+
+Sculley, D. (2010). Web-scale k-means clustering. *Proceedings of the 19th International Conference on World Wide Web*, 1177--1178.
+
+Strehl, A. & Ghosh, J. (2002). Cluster ensembles -- a knowledge reuse framework for combining multiple partitions. *Journal of Machine Learning Research*, 3, 583--617.
+
+Van Gassen, S., Callebaut, B., Van Helden, M.J., et al. (2015). FlowSOM: Using self-organizing maps for visualization and interpretation of cytometry data. *Cytometry Part A*, 87(7), 636--645.
+
+Weber, L.M. & Robinson, M.D. (2016). Comparison of clustering methods for high-dimensional single-cell flow and mass cytometry data. *Cytometry Part A*, 89(12), 1084--1096.
+
+Wilkerson, M.D. & Hayes, D.N. (2010). ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking. *Bioinformatics*, 26(12), 1572--1573.