arXiv CSopen access

Cluster LOCO: Feature Importance For Interpreting Clusters

Cluster LOCO: Feature Importance for Interpreting Clusters

Cluster LOCO: Feature Importance For Interpreting Clusters Claire He

[email protected]

Department of Statistics Columbia University New York, NY 10027, USA

arXiv:2606.14592v1 [stat.ML] 12 Jun 2026

Genevera Allen

[email protected]

Department of Statistics Columbia University New York, NY 10027, USA

Abstract Clustering is widely used for exploratory analysis and scientific discovery, driving insights from market segmentation to biological data analysis, but its outputs can be difficult to interpret, audit, and reproduce as modern datasets become increasingly large and complex. Reliable use of clustering requires understanding which features drive the discovered structure, yet feature-level explanations for clustering remain scarce compared with methods in supervised learning. Furthermore, existing clustering feature importance scores are often tied to specific algorithms and data assumptions. To address these challenges, we propose Cluster LOCO (Leave-One-Covariate-Out), a family of model-agnostic feature importance scores for clustering. Cluster LOCO is built on feature occlusion and clustering generalizability, defined as whether cluster labels learned on one subset of the data can be accurately predicted on held-out samples. For any chosen clustering algorithm, Cluster LOCO quantifies a feature’s importance by measuring how much its removal degrades generalizability. We first introduce Cluster LOCO-Split, which relies on data splitting, and then extend it to Cluster LOCO-MP, a minipatch ensemble-based version designed for large-scale data. Across synthetic simulations and an application to cell-type discovery in single-cell transcriptomics, we show that Cluster LOCO more reliably recovers informative features than existing clustering feature importance methods.

1 Introduction A fundamental task in unsupervised learning, clustering is used across disciplines ranging from the social sciences, astrophysics to biology [17; 41; 66] to draw insights from data by forming groups or partitions. Yet clustering is not defined by one canonical objective. As emphasized by Luxburg et al. [38], clustering has several use cases: it may be used for preprocessing, to organize, compress, or denoise data; for exploration, to reveal unknown structure and generate hypotheses; or for confirmation, to validate hypothesized groupings or support scientific discoveries. Because clustering intervenes upstream in the data science life-cycle [70], its outputs influence downstream analysis, modeling, and scientific interpretation. It is therefore consequential that clustering-driven conclusions be reliable, trustworthy, and reproducible. At the same time, reliability is difficult to assess for clustering solutions because of the underlying assumptions and choices made by clustering models: for example, K-means re1

quires specifying the number of clusters, and the algorithm will return exactly that many groups whether or not such structure is meaningful in the data [2]. This challenge is amplified in modern datasets, where complex nonlinear patterns, interactions, and high dimensionality can limit the effectiveness of classical methods [26]. In response, practitioners increasingly rely on deep clustering models [42; 32] or heavily feature-engineered workflows [8; 22; 21; 53; 63]. While these approaches can improve the detection of complex structure (e.g. nonlinearity) or domain-dependent specificity (e.g. zero-inflated data), they can also make the resulting clusters harder to understand and audit. For instance, in genomics, clustering has enabled meaningful discoveries of cell types and markers [60] while, on the other hand, computational studies centered on clustering have raised persistent concerns about reproducibility [16]. For clustering to support rigorous discovery, we therefore need tools that clarify why a clustering solution arises. In particular, feature-level explanations can help identify which features the clustering solution relies on and promote trust [15]. We address this need by bringing the perspective of interpretable machine learning (IML) to clustering, establishing a useful notion of feature importance for clustering solutions. In many clustering applications, features are themselves meaningful and interpretable: they are the genes in transcriptomic data, words in text data, behavioral or measured attributes in the social sciences and often the object of downstream analysis. For example, in single-cell genomics, practitioners commonly interpret clusters by identifying “marker genes” i.e. genes that differ across the discovered clustering groups [60], through differential expression analysis, where genes are tested across clusters [25]. While useful for annotation, this workflow can raise post-selection inference, or “double dipping” concerns because the same data are used both to define the clusters and to assess the features that distinguish them, a practice known to inflate the false discovery rate of significant genes [71; 7; 52]. Moreover, features that differ across clusters are not necessarily the features that produced the clustering solution. We therefore focus on feature importance for clustering, which asks directly how much each feature contributes to the clustering structure itself. This IML perspective is well established in supervised learning, where feature importance methods are widely used to explain model predictions [43]. In clustering, however, featurelevel interpretability remains comparatively underdeveloped. Existing work has largely followed two directions: intrinsically interpretable clustering algorithms, built on interpretable supervised methods such as decision trees with a modified objective for clustering [20], and post-hoc explanations tailored to specific algorithms, mainly K-means [44; 24]. While useful, these approaches can be difficult to scale to large datasets, might make particular data assumptions or constrain the practitioner to specific clustering models. On the other hand, feature selection in clustering enables handling large scale data and has been well studied with models leveraging sparsity via regularization [62; 61], filters or wrapper methods [65; 6; 49; 1]. While related, feature selection and feature importance answer different questions. Feature selection asks which features should be used to construct a clustering solution, often by optimizing a sparsity or clustering quality criterion. In contrast, feature importance asks, after a clustering solution has been obtained, which features contributed most to it. We instead seek a post-hoc interpretation of a chosen clustering model, separat2

Cluster LOCO: Feature Importance for Interpreting Clusters

ing the explanation of a clustering solution from the decisions used to produce it. To formulate such a post-hoc interpretation, we turn to feature importance methods from supervised interpretable machine learning. In supervised learning, feature-level interpretability includes a wide range of model-specific and model-agnostic techniques. For model-agnostic techniques, we find three main types of feature importance metrics: feature permutation introduced via model-class reliance [12; 4] inspects the model’s "reliance" to each feature via its error; Shapley values [51; 36] and extensions [54; 40; 59], a popular metric based in game-theoretical axioms which distributes feature value across features fairly; and feature ablation or occlusion which explains the change in prediction when a feature is removed [29; 48; 59]. In this paper, we propose a model-agnostic feature importance framework for clustering based on feature occlusion. Our approach is motivated by the Leave-One-Covariate-Out framework [29]: remove a feature and measure how much the clustering solution changes. If removing a feature substantially alters the clustering, then the feature is important for the discovered structure; if the clustering remains largely unchanged, then the feature is less important. This definition is straightforward to interpret, post hoc as it is applied after a clustering method has already been chosen, and is therefore flexible to complex clustering workflows. Our contributions are as follows: first, we introduce Cluster LOCO, a family of model-agnostic feature importance metrics for clustering based on feature occlusion via Cluster-LOCO-Split. Second, we develop Cluster LOCO-MP, a scalable minipatch-based extension designed for high-dimensional data. Third, we validate our approach in simulated settings, including low-dimensional examples with complex nonlinear structure and highdimensional regimes, and compare it against existing feature importance methods. Lastly, we demonstrate our framework in a single-cell transcriptomics application where important features for clustering are biologically meaningful.

2 Feature importance scores for clustering via LOCO 2.1 Background for Cluster LOCO Before introducing our Cluster LOCO family of metrics, we review two existing notions that inspired our metric: generalizability and LOCO. Generalizability in clustering, also sometimes coined as predictability, was first proposed in the context of model selection by Lange et al. [28] and Tibshirani and Walther [56] as an alternative to stability-based validation techniques relying on sampling or bootstrapping non disjoint subsets of data [3; 31]. Instead, the stability of cluster solution is captured by a transfer predictor to measure how generalizable the cluster solution on one disjoint subset is to another. We apply this idea to the ML paradigm of data splitting, where a portion of the data is carved out as a training set, another held-out for calibration i.e. X = Xtr ⊔ Xcal . Then the measure of generalizability of the clustering solution is defined as the error (or dissimilarity) between the clustering solution on the calibration set (zcal cluster labels for Xcal ) with its transferred labels via a transfer classifier fˆtr trained on the 3

training data (Xtr with cluster labels ztr ). generalizability error = Error zcal , fˆtr (Xcal ; (Xtr , ztr )

(1)

This quantifies the generalizability of clustering solutions: a small value for the error (or dissimilarity) between clustering solutions and predictions on the calibration set means the clustering solution learned from the train set is very generalizable to the held-out calibration set. Example of error measures include the mean squared error (used observation-wise in [56]’s prediction strength index), or the Hamming distance (which corresponds to the misclassification risk of the calibration set) in [28]. However, for clustering which involves multi-class labels, a more popular similarity metric for clustering label comparison is the adjusted-Rand index (ARI). In fact, the aforementioned error measures require a label alignment step, a linear assignment problem that can be solved using the Hungarian matching algorithm [27] while ARI is permutation invariant. Alternatively, we also use the multi-class hinge loss as point-wise error measure for Cluster LOCO when using a soft classifier. These scores empirically yield similar normalized feature importance (see Appendix Figure 5 for a comparison of error measures on a simple example). Leave-One-Covariate-Out (LOCO) on the other hand is an extensively studied quantity borrowed from supervised learning [29; 59; 14; 33] which quantifies the change in prediction error when removing a feature. In its LOCO-Split form, the LOCO error for feature j given test point (X, y) is obtained by the difference in prediction error on the test set in the −j absence of feature j (i.e. fˆtr fit on without-j training data (Xtr,−j , ytr )) and the prediction error on the test set given the full feature set for training (i.e. fˆtr trained on (Xtr , ytr )) as shown in equation 2.

−j LOCOj (X, y) = Error y, fˆtr (X−j ; (Xtr,−j , ytr ) − Error y, fˆtr (X; (Xtr , ytr )

(2)

Essentially, LOCO quantifies how much the model performance drops when retraining the model after excluding feature j to determine whether it was an important feature: if the performance degrades, then feature j is important whereas if the performance remains unchanged, feature j isn’t important. 2.2 Cluster LOCO-Split: a generalizability feature importance score The LOCO objective lends itself naturally to a clustering extension: one may remove a feature, re-cluster the data, measure the resulting change in cluster solution using internal validity indices or stability scores for example. However, popular validity indices such as the silhouette score [50; 23] often rely on geometric assumptions (i.e. favoring compact or well-separated clusters). On the other hand, most stability criteria [3; 37] measure the robustness of the clustering solution under resampling and the LOCO interpretation of the effect of feature occlusion is less meaningful as scores reflect mostly the model’s sensitivity. Therefore, we choose to use generalizability as a meaningful quantity for feature importance via LOCO: providing a model-agnostic and assumption-free feature importance. We give an illustrative example of both silhouette-based LOCO and stability-based LOCO scores’ 4

Cluster LOCO: Feature Importance for Interpreting Clusters

limitations in Section 3 compared to our generalizability-based Cluster LOCO-Split score. In order to capture the change in generalizability of clustering solutions due to feature contribution, our score (in eq. 3) quantifies the change in generalizability error when removing a feature. Essentially, when an important feature is removed for the cluster algorithm, the clustering solution is expected to be less generalizable and the generalizability error without the feature increases. Cluster LOCO-Split. For X ∈ RN ×M split into Xtr and Xcal , Cθ a clustering algorithm. Cluster the training set and get the cluster labels ztr = Cθ (Xtr ), cluster the calibration set and get cluster labels zcal = Cθ (Xcal ):

ˆ j (X) = Error zcal , fˆ−j (Xcal,−j ; Xtr,−j , ztr ) − Error zcal , fˆtr (Xcal ; (Xtr , ztr ) ∆ tr | {z } | {z } generalizability error without j

(3)

generalizability error

We summarize our Cluster LOCO-Split algorithm in Algorithm 1, where the metric of dissimilarity used from the stability literature is usually the negative ARI. We start by splitting the data into training and calibration i.e. Xtr = {Xi }i∈Itr and Xcal = {Xi }i∈Ical respectively, and for cluster algorithm Cθ (where θ denotes all hyperparameters for the clustering algorithm), we separately obtain the cluster labels on each split ztr = Cθ (Xtr ), zcal = Cθ (Xcal ). We set {Xtr , ztr } to be the training set for the generalizability predictor while {Xcal , zcal } is used as unseen held-out set. Cluster LOCO-Split requires two generalizability classifiers, a base generalizability classifier fˆtr trained on the full training set (Xtr , ztr ) and its without −j feature-j counterpart fˆtr trained on the without feature-j training set (Xtr,−j , ztr ). We derive the Cluster LOCO-Split score as the difference between the generalizability error leaving out feature j evaluated on the calibration set (Xcal,−j , zcal ) and the generalizability error with full features evaluated on the full feature calibration set (Xcal , zcal ). Algorithm 1 Cluster LOCO-Split Require: Unlabeled X ∈ RN ×M 1: Split data into training and calibration sets: Xtr = {Xi }i∈Itr and Xcal = {Xi }i∈Ical . 2: Cluster each split: (a) Cluster the training data: {Xi }i∈Itr to obtain cluster labels {zi }i∈Itr . (b) Cluster the calibration data: {Xi }i∈Ical to obtain cluster labels {zi }i∈Ical . 3: Fit generalizability classifiers (a) Fit fˆ to the cluster-labeled training data {(Xi , zi )}i∈Itr . (b) Fit fˆ−j to cluster-labeled without j training data {(Xi,−j , zi )}i∈Itr . 4: Compute & return Cluster LOCO-Split

ˆ split := Error zcal , fˆ−j (Xcal,−j |(Xtr,−j , ztr )) − Error zcal , fˆtr (Xcal |(Xtr , ztr )) ∆ tr j ˆ split }M Ensure: {∆ j=1 j

5

2.3 Scaling Cluster LOCO: a fast procedure for high-dimensional data Cluster LOCO-Split inherits several limitations from data splitting and feature occlusion. First, splitting the data reduces the sample size available both for fitting the clustering model and for estimating the generalizability error. This can reduce accuracy and make the score sensitive to the specific train-calibration split. The issue is especially pronounced when clusters are unbalanced: if either split contains few observations from a cluster, the resulting generalizability estimate can become unstable. Second, because Cluster LOCO-Split removes one feature at a time, it can underestimate the importance of correlated features: in the presence of two correlated features, removing one may have little effect because the other remains in the active feature set; as a result, both features may receive artificially low importance scores. Multi-fold splitting can mitigate some of the instability caused by data splitting, but it increases the computational cost substantially and may be impractical for large datasets as one would have to refit F × (M + 1) models, where F the number of folds. In this section, we introduce Cluster LOCO-MP, a minipatch-based extension designed to address these limitations. By leveraging ensembles of small random subsets of observations and features, Cluster LOCO-MP provides a fast and more stable approximation of Cluster LOCO scores, improving scalability in large datasets.

2.3.1 Cluster LOCO-MP: Cluster LOCO with minipatch ensembles To scale up our Cluster LOCO feature importance score for large-scale data, we introduce minipatches. Minipatches are tiny subsamples of the data that enable fast and extremely parallel model fitting via ensembling. Instead of traditional subsampling observations into batches [4; 34], minipatches require simultaneous subsampling of n << N observations and m << M features, and yield a natural structure of the data in in- and out-of-patch features and in- and out-of-patch observations. Minipatches have been shown to have computational advantages in the supervised learning setting [67; 57; 14], benefit from implicit ridge-like regularization useful in correlated settings [68] and have been used successfully for clustering in consensus clustering [13]. By drawing on those properties of minipatches, our Cluster LOCO-MP score is fast with large datasets, and no longer sensitive to correlated features. Algorithm 2 presents the full Cluster LOCO-MP procedure. The algorithm follows the same principle as Cluster LOCO-Split, but replaces a single train-calibration split with an ensemble of minipatches in a leave-one-out framework exploiting the natural in and outof-patch structure. Given a base cluster algorithm and classification algorithm, we require to fix ahead a large number of minipatches B. Then our algorithm can be summarized in four steps: first, performing minipatch clustering and generalizability training; second and third, obtaining the ensemble LOO and LOCO-LOO predictions; and finally computing the scores. Since we cluster each minipatch and then train a generalizability classifier to predict the resulting cluster labels, they are arbitrary across minipatches and require to be aligned to a common labeling. We outline in our algorithm an alignment step using a reference clustering but when computing such a reference clustering is costly, one can alternatively use the pairwise overlap across minipatches to get approximate alignment. The aligned predictors are then aggregated out-of-training: the LOO ensemble predictor Ĥ averages 6

Cluster LOCO: Feature Importance for Interpreting Clusters

soft-predictions (or uses max-vote with hard-predictions) over minipatches that exclude the target observation, whereas the LOCO-LOO predictor Ĥ −j further restricts this average to minipatches that also exclude feature j. The Cluster LOCO-MP score is then defined as the change in prediction error between the ensemble Ĥ −j and Ĥ. Algorithm 2 Cluster LOCO-MP Require: Unlabeled X1 , . . . , XN ∈ RM to cluster. 1: Perform minipatch consensus clustering and minipatch training: (a) Randomly sample B minipatches of size n × m indexed (Ib , Fb )B b=1 . (b) Obtain a reference clustering z0 = Cθ (X) for alignment. (c) Cluster each minipatch b to get the minipatch cluster label z̃Ib = Cθ (XIb ,Fb ). Align z̃Ib with z0,Ib via Hungarian matching to get zIb . (d) Calculate consensus clustering z ∗ where for i = 1, . . . , N , X zi∗ = arg max 1(i ∈ Ib )1(zi = k) k∈[K]

b∈[B]

(e) Train the generalizability classifiers on each minipatch b: fˆb on (XIb ,Fb , zIb ). 2: Construct Cluster LOO ensemble prediction Ĥ for observation i: X Ĥ(Xi,: ) := fˆb (Xi,Fb |(XIb ,Fb , zIb )) {b∈[B]:i∈I / b}

3: Construct Cluster LOCO-LOO ensemble prediction Ĥ −j for features j = 1, . . . , M and

observation i:

Ĥ −j (Xi,−j ) :=

X

fˆb (Xi,Fb |(XIb ,Fb , zIb ))

{b∈[B]:i∈I / b ,j ∈F / b}

4: Compute Cluster LOCO-MP scores for feature j = 1, . . . , M :

ˆ j := Error(z ∗ , Ĥ −j (X:,−j )) − Error(z ∗ , Ĥ(X)) ∆ ˆ j }M Ensure: {∆ j=1 The computational advantages of our method come primarily from the minipatch construction. Each minipatch contains only a small subset of observations and features, reducing the cost of fitting computationally expensive clustering procedures on large data. Moreover, because minipatches are generated and fitted independently, the training step is naturally parallelizable. For example, spectral clustering has time complexity O(N 3 ), which becomes O(Bn3 ) under minipatch clustering, where n ≪ N is the number of observations in each minipatch and B is the number of minipatches. The ensemble scores can also be computed in parallel across features, making the implementation modular and adaptable to the available computational resources. Beyond computational efficiency, feature subsampling in the minipatch ensemble helps evaluate the contributions of correlated features. Because different 7

minipatches contain different subsets of features, some include one feature of a correlated group but not another, while others exclude the group altogether. The resulting ensemble prediction error therefore reflects how much predictive information is lost when a feature is unavailable across the ensemble, even when it belongs to a set of correlated features.

2.3.2 Cluster LOCO-RAMPART In high-dimensional applications, practitioners are often primarily interested in the top-k ranked features driving a clustering solution rather than the full set of important features. This is particularly relevant in applications such as genomics, where many features may be irrelevant or redundant (e.g. house-keeping genes) and interpretation typically focuses on a small set of candidate markers. In such cases, computing Cluster LOCO-MP scores for every feature can be both computationally expensive and statistically inefficient. We therefore combine Cluster LOCO with the RAMPART framework of Chen et al. [5], a data-adaptive procedure for top-k feature ranking. In RAMPART Cluster LOCO, Cluster LOCO-MP is used as the round-level feature scoring oracle: at round t, Cluster LOCO-MP produces imˆ (t) : j ∈ St } on the active feature set St . Features are then ranked from portance scores {∆ j these scores, the bottom half is discarded, and the procedure repeats until only the topk candidates remain. This yields an adaptive feature importance procedure: preliminary Cluster LOCO-MP estimates are used to screen and prioritize candidate features, and additional computation is concentrated on the most promising features retained through multiple rounds of fixed-batch sequential halving. Chen et al. [5] showed that while theoretically the order of the number of minipatches required with this procedure is roughly the same as for the minipatch estimates, empirically, in their work and also in our sparse high-dimensional simulations in section 3, this adaptive procedure improves the identification of important features while reducing computational cost substantially. The underlying interpretation remains unchanged: features are important when their removal increases the generalizability error of the clustering solution. The full Cluster LOCO-RAMPART algorithm is provided in the Appendix (Algorithm 3).

2.4 Extensions and practicalities 2.4.1 Cluster-specific LOCO In many applications, cluster-level feature interpretability is more informative than global feature interpretability: for instance, in genomics, this corresponds to understanding which genes are important for a particular cluster that might identify with a specific cell type. We defined earlier the Cluster LOCO importance metric in the most general sense, compatible with dissimilarity metrics as well as point-wise error scores. In the case of point-wise error scores like multi-class hinge loss, we can extend the family of Cluster LOCO global importance scores to cluster-level importance scores. Cluster LOCO-Split admits a natural cluster-specific interpretation: instead of aggregating over all data points, one can aggregate over each cluster and get a cluster-specific important score (see eq. 4). For feature j and 8

Cluster LOCO: Feature Importance for Interpreting Clusters

cluster k, let Nkcal = ˆ split = ∆ j,k

P

i∈Ical 1(zi = k):

X 1(zi = k)  i∈Ical

Nkcal

−j Error zi , fˆtr (Xi,−j ; (Xtr,−j , ztr )) −Error zi , fˆtr (Xi,: ; (Xtr , ztr ))



(4) Similarly, we can derive the Cluster-specific LOCO-MP estimates for any point-wise error metric, where for feature j and cluster k, ˆ MP = ∆ j,k

N X i=1

1(zi∗ = k) PN

∗ i=1 1(zi = k)

 

Error zi∗ , Ĥ −j (Xi,−j ) − Error zi∗ , Ĥ(Xi,: )

(5)

This property makes our score highly interpretable at both a local and global scale which is not the case for most scores in the literature: only the SHAP-derived quantities and LRP admit sample-level to global-level interpretation. For clustering however it is especially important to be able to assess the model at the cluster level, especially when those clusters are used for scientific analysis in downstream tasks. We present such a case in section 4. 2.4.2 Choice of hyperparameters Cluster LOCO is a family of post-hoc interpretability metrics and is not intended to replace careful selection and validation of the clustering model itself. The choice of clustering algorithm, distance or similarity measure, number of clusters, and other model-specific hyperparameters should be made before using Cluster LOCO metrics, drawing on appropriate domain knowledge and clustering validation procedures (see Yin and Hamerly [69]; Luxburg [37]; Handl et al. [18]; Ullmann et al. [58]; Gan et al. [15]. Once a clustering solution has been selected, Cluster LOCO metrics can then be used to interpret which features contribute most to that solution. For Cluster LOCO-Split and Cluster LOCO-MP, the main tuning choices are the train/calibration split ratios and minipatch size ratios respectively. For the former, both splits must contain enough observations to obtain meaningful cluster labels and stable generalizability estimates. We recommend shuffling the data before splitting and using a balanced split, such as a 50% train and 50% calibration split, as a default. This choice helps avoid pathological splits in which rare or unbalanced clusters are poorly represented in either the training or calibration set. For Cluster LOCO-MP, the main hyperparameters are the minipatch size ratios. The minipatch size controls the tradeoff between computational efficiency and the accuracy of the score estimates: smaller minipatches are faster but may yield noisier clustering solutions, while larger minipatches better approximate the full-data clustering problem at greater computational cost. The minipatch size ratios can be tuned in general as discussed in Gan et al. [14] and tend to be robust when selecting the appropriate n m ratios (rN = N , rM = M ) in the range of 20%-50% for both observations and features [13].

3 Simulations In this section, we evaluate our proposed Cluster LOCO-Split and Cluster LOCO-MP/LOCORAMPART methods against existing feature importance approaches in simulation settings. We compare our metrics against five baseline methods that provide feature importance for clustering: a prototype-based feature importance (PBFI) and Fuzzy C-Means Shapley (FCM 9

SHAP) from Nápoles et al. [44], permutation feature importance (PFI) [46], an extension of layer-wise relevance propagation (LRP) for KMeans clustering proposed by Kauffmann et al. [24]; and IMPACC, a minipatch-based adaptive consensus clustering algorithm proposed by Gan and Allen [13]. Since no public code base was available for PBFI or FCM SHAP, we re-implemented both methods in Python. We also extended the KMeans-based PFI implementation, originally proposed in R, to arbitrary scikit-learn-style clustering algorithms, and similarly adapted IMPACC, originally in R as well, to be compatible with these algorithms. 3.1 Illustrative simulation for Cluster LOCO-Split To illustrate the limitations of existing methods and how Cluster LOCO-Split is able to obtain better feature importance scores in even simple examples, we construct a small but challenging synthetic example for clustering. This synthetic clustering example consists of 3 signal features created from interlaced half circles with varying levels of overlap and noise shown in Figure 1a; details on the data-generating process are provided in the Appendix. We augment those 3 signal features with 2 pure noise features sampled from independent uniform distributions. The difficulty of this example stems from the nonlinear separability of the clusters, which challenges methods that favor convex cluster geometries. Spectral Clustering recovers the true groups as indicated by an adjusted Rand index of 1, whereas methods that impose convex or linearly separable clusters, such as KMeans with an adjusted Rand index of 0.40, produce incorrect assignments as shown in Figure 1b. Consequently, feature importance methods using KMeans, including FCM-SHAP (where we used exact Shapley), LRP, PBFI, and the original permutation-based score PFI, explain a misspecified clustering model rather than the clustering structure of interest. As shown in Figure 1b, these methods assign high importance to noise features in this setting. We also investigated several LOCO-style scores based on different clustering validation criteria: we show two such variants constructed from popular scores, silhouette LOCO is based on the silhouette score change when removing one feature, stability LOCO computes the change in negative ARI under the model-explorer stability objective, using the framework of [3]. We find that our generalizability-based score, Cluster LOCO-Split, provides the most informative notion of feature importance in this simple yet challenging setting where it uniquely separates signal from noise and correctly orders the signal features by importance.

3.2 Cluster LOCO for large data via minipatch ensembles We next applied Cluster LOCO-MP for large data: generating three base clustering simulations with K = 7 clusters, N = 3500 observations (500 observations per cluster), with fixed signal features p∗ = 10 and increasing noise features pnoise ∈ {10, 40, 190, 490, 990} covering different clustering difficulties. We propose two mixtures models: a Gaussian mixture with onion covariance structure (as proposed by Qiu and Joe [47]) to introduce correlation in the features and a Gamma mixture model that generates data with heavier tails. For a more complex scenario of clustering, we generate clusters inspired by the two-dimensional toy examples of moons and concentric circles [10]: we sample our data in a two-dimensional 10

Cluster LOCO: Feature Importance for Interpreting Clusters

a.

X2

X3

X3

X1 b.

X5

X2

X1

Cluster LOCO Split

Spectral clustering

Stability LOCO (Spectral Clustering)

(Spectral Clustering)

X4 Silhouette LOCO (Spectral Clustering) Method Spectral KMeans Feature type Signal Noise

X2

LOCO scores

X1 X2 X3 X4

ARIKMeans = 0.40 ARISpectral = 1.00

X5 0.00

X2

Benchmark methods

X1 KMeans

X1

0.05 Permutation FI (KMeans)

0.10

0.2 0.1 0.0 0.1 Shapley (permutation) (Fuzzy c-Means)

0.00

0.02 0.04 LRP (neuralized KMeans)

Prototype FI (KMeans)

X1 X2 X3 X4 X5 0.0

0.2

0.4 0.00

0.02

0.04

0.06 0.00

0.05

0.10

0.0

0.1

0.2

0.3

Figure 1: Cluster LOCO is a reliable feature importance score for complex clustering problems. a. Simulated dataset presents nonlinearly separable clusters, first three features are signal features, ordered by their importance of contribution. Removing X1 leads to a harder clustering problem in the remaining signal space with overlapping nonlinearly separable clusters, removing X2 leaves a moderately hard problem while removing X3 leads to the easier problem hence the ordered importance of each feature. The last two features are pure noise features sampled from a uniform distribution. b. Feature importance scores of possible LOCO-style metrics and existing IML methods for clustering in the literature. For model-agnostic methods we use Spectral Clustering that yields the highest adjusted Rand index (ARI) on this problem while we report the method used for model-specific algorithms. Our Cluster LOCO importance score is the only metric able to recover the correct signal features with the correct ordering. embedding with controlled geometry and project up using an orthogonal projection to get any higher-dimensional extension. Each of those three base simulations has two difficulty levels created by varying cluster separation or signal-to-noise ratio into easier or harder clustering problems as shown in the PC-space scatter plot of the generated data in Figure 2 row 1 and 3. We report the top-10 hits (or precision at 10 features) of the true 10 signal features 1 for each method i.e. 10 |{top 10 features by method} ∩ {true 10-feature set}|. Our Clus11

B1.

C1.

Gamma Mixture (easier)

Interlaced Moons/Circles (easier)

PC 2

PC 2

Gaussian Mixture (easier)

PC 2

A1.

PC 1

PC 1

PC 1 C2.

1.0

1.0

0.8

0.8

0.8

0.6 0.4

0.6

Top-k Hits

Top-k Hits

B2.

1.0

Top-k Hits

A2.

0.4 0.2

0.2

0.6 0.4 0.2

0.0 20

200

Number of Features p

500

1000

20

B3.

50

200

Number of Features p

500

20

1000

PC 1

1.0

0.8

0.8

0.8

0.6

0.6

Top-k Hits

C4.

1.0

Top-k Hits

B4.

1.0

Top-k Hits

A4.

0.4

0.4

0.0 50

200

Number of Features p

Cluster LOCO RAMPART

500

1000

Cluster LOCO-MP

0.4 0.2

0.2

0.2

500 1000

Interlaced Moons/Circles (harder)

PC 1

0.6

200

Number of Features p

PC 2

PC 1

20

50

C3.

Gamma Mixture (harder)

PC 2

Gaussian Mixture (harder)

PC 2

A3.

50

20

50

200

Number of Features p

PBFI (KMeans)

500

LRP (KMeans)

1000

IMPACC

0.0

20

50

200

Number of Features p

Permutation

500 1000

Fuzzy C-Means SHAP

Figure 2: Simulation results are obtained across three settings: A. Gaussian mixture, B. Gamma mixture, and C. interlaced moons/circles. Model-agnostic algorithms used are KMeans for Gaussian mixtures, Gamma mixture EM for Gamma mixture, Spectral Clustering for interlaced moons/circles. For each setting, rows 1 and 3 show the data colored by true labels in PC space, while rows 2 and 4 report top-10 hits for feature importance with signal features p∗ = 10 and noise features pnoise ∈ {10, 40, 190, 490, 990}. 12 Model-specific feature importance scores’ models are reported in the legend. Top-k hits were averaged over 100 replicates. Cluster LOCO-MP and RAMPART outperform existing methods on most tasks.

Cluster LOCO: Feature Importance for Interpreting Clusters

ter LOCO-MP and Cluster LOCO RAMPART methods were run with fixed B = 5000 and Brampart = 1000 and minipatch ratio sizes αN = 0.2, αM = 0.2, while IMPACC was run with the default hyperparameters. Because FCM SHAP is extremely computationally expensive to run for exact Shapley value, we used in this example the SHAP permutation approximation of Shapley values with 20 iterations. For each simulation setting, model-agnostic methods were evaluated with the best performing algorithm among typical clustering algorithm of the scikit-learn environment: for Gaussian mixtures, we used KMeans, for the interlaced moons/circles, we used Spectral Clustering. For Gamma mixture, we implemented an EM Gamma mixture model that was best suited for clustering this type of data. In Figure 2 panels A2 and A4, the performance of the different feature importance scores in the Gaussian mixture model case is reported via the top-10 hits: PBFI and LRP which are KMeans-based recover the correct features as expected, our Cluster LOCO-RAMPART and Cluster LOCO-MP scores also recover the correct top-10 features for p ≤ 500, with a slight drop in performance at p = 1000 which could be solved by budgeting ahead more minipatches. In the easy case, IMPACC performs well, but improves especially when the number of noise features grows as it is designed for sparse features. FCM SHAP and PFI perform the worst overall in the easy and hard case of the Gaussian mixture. For FCM SHAP, the steep drop off of performance is due to the gap between exact Shapley and the approximation with fixed number of iterations while the number of features grow. In the two other simulation examples in Figure 2 columns B and C, our Cluster LOCO family scores outperforms all comparative scores the most definitely: in the Gamma mixture case, IMPACC achieves similar or better performance only when the features are extremely sparse (5% or less of the features are signal features), while in the easy case, both Cluster LOCO RAMPART and Cluster LOCO-MP have robust recovery of the top-10 features. We see similar results with the interlaced moons/circles (Figure 2 panels C2 and C4), Cluster LOCO-MP and Cluster LOCO RAMPART being extremely consistent and effective at recovering the true signal features. We note that in the easy case, LRP is also able to recover effectively the true signal when there is stronger separation in the data and nonlinearity does not affect KMeans as much when it was initialized well. However this performance drops when the clusters have more overlap in their supports (see Figure 2 panels C3 and C4). Overall, Cluster LOCO methods for large data perform as well as or outperform existing methods for those diverse clustering simulations, even typically "hard" clustering tasks, and we provide further simulation results in the Additional figures.

4 An application in single-cell transcriptomics to immune cells We apply Cluster LOCO-MP and Cluster LOCO-RAMPART for scientific discovery in a real-world clustering application. In single-cell transcriptomics, clustering is routinely used to discover putative cell types or states, after which differentially expressed genes are identified post-hoc and interpreted as marker genes for the discovered groups. This workflow is embedded in popular pipelines such as Scanpy and Seurat [63; 19], and sometimes considered the gold standard for discovering cell-type identities [35]. However, this two-step procedure raises important statistical concerns: the same data are first used to define clusters and then reused to test for genes that distinguish those clusters, a form of double dipping or data 13

snooping known to inflate false discovery rates [39]. Moreover, this sort of cell-type clustering for gene marker annotation is not a unified framework: the clustering algorithm choice can substantially change the downstream discoveries, and those new marker genes are implicitly dependent on the clustering solution, which can contribute to a lack of reproducibility of results [16]. Cluster LOCO addresses this gap by directly quantifying the contribution of each feature to clustering generalizability, providing a clear unified framework for feature importance in clustering rather than relying on separate downstream testing steps. 4.1 Peripheral Blood Mononuclear Cells dataset for immune cell types discovery We study the Peripheral Blood Mononuclear Cells (PBMC) dataset published by Zheng et al. [72] that was annotated using purified transcriptomics populations. Our goal is to cluster and understand the cluster-wise marker genes discovered using traditional pipelines of gene annotations (known to suffer from double dipping) and our proposal of feature importance. We use the processed data available via Scanpy [63] that consists of 765 genes and 700 single cells, obtained after normalizing and scaling the data as reported on the 10X Genomics repository. The reported labels were obtained by classifying the single-cell data with 11 purified sub-populations of PBMC reference profiles and are referred to as bulk labels or purified labels. It is important here that we chose data that was not generated via clustering but via a correlation approach for external validation of our results. To further illustrate our results, we compiled known marker annotations of human PBMCs from three sources: the Azimuth atlas [53], as well as Ding et al. [9] and Oelen et al. [45] supplementary data on known markers of human PBMC cell types. When pooling known markers, we note that the granularity of cell typing was not always consistent, and this leads to some markers ambiguously marking multiple cell types. 4.2 Globally important genes for clustering In this data application, global feature importance scores identify the genes most influential for the overall clustering solution. A reliable and interpretable model should reflect the biological reality, therefore we expect highest ranked features to align with some of our known marker genes defining cellular identity in human PBMC. We report in Figure 3 the results for six clustering feature importance metrics: Cluster LOCO-MP, Cluster LOCO-RAMPART with top-100 genes, IMPACC (consensus clustering with minipatches importance), LRP (layer-wise relevant propagation score), PBFI (prototype-based feature importance) and PFI (permutation feature importance). Following the best practices for model selection in clustering [2; 64], we chose hierarchical clustering with Ward linkage with K = 10 clusters to correspond with the 10 known purified labels as base clustering algorithm for our model-agnostic feature importance scores. We show in Figure 3a the clustering solutions in PC-space aligned with the reported labels for both hierarchical clustering and KMeans clustering that PBFI and PFI use. Cluster LOCO-MP was run with B = 5000 minipatches and minipatch ratio sizes αN = 0.42 for observations, αM = 0.26 for features. We aligned the obtained labels for each methods with the reported purified labels using the Hungarian algorithm.

14

Cluster LOCO: Feature Importance for Interpreting Clusters

To evaluate the different clustering feature importance methods, we first analyze the top 10 most globally important features in Figure 3b, where the highest-ranked features correspond to the most important genes driving the clustering structure. We first note that KMeansbased importance methods PBFI and PFI fail to find any known marker genes among the 10 most important genes. Among the remaining four methods, the CD14+ monocyte marker FTL is consistently important for Cluster LOCO-MP, Cluster LOCO-RAMPART, LRP and IMPACC, matching the original empirical findings in Zheng et al. [72]. Overall, Cluster LOCO-MP recovers a larger number of reference cell type-specific marker genes than the competing methods. Cluster LOCO-MP also demonstrates the strongest alignment with the biological ground truth, assigning high importance scores primarily to known marker genes from dendritic cells and monocytes, which correspond to the two most well-separated clusters in the PC space. On the other hand, LRP that uses a neuralized KMeans and therefore is more flexible than traditional KMeans recovers similar reference markers as Cluster LOCO-MP. In particular, Cluster LOCO-MP and LRP share important genes like AIF1 and LST1 that do not belong in our reference marker set but whose expression patterns have been documented in the monocyte and myeloid cell types in the PBMC literature [55; 30; 11]. However, LRP also presents PSAP a gene with lower immune-lineage specificity [55], and therefore with a less interpretable profile for cell typing. We note that IMPACC’s topranked features include a mixture of CD34+ hematopoietic stem cell markers and monocyte markers. However, half of its top 10 most important genes lack documented cell-type specificity (e.g., IGLL1 ). These results underscore that in order to interpret important genes, the clustering solution needs to match a relevant biological truth. In particular, selecting hierarchical clustering as opposed to KMeans yields clustering solution that match a closer biological truth and therefore the solutions are necessarily more interpretable. 4.3 Cluster-wise important genes and marker selection Global feature importance scores identify the genes most influential for the overall clustering solution, but cluster-level scores are needed to assess which genes characterize individual cell types. We compare here Cluster LOCO-MP at cluster-level with LRP also aggregated at cluster-level – other importance scores investigated earlier are solely global feature importance scores. In parallel, we implement the clustering + differential gene expression workflow commonly used in single-cell analysis and described in Scanpy tutorials [63]: data is clustered using Louvain clustering, identifying 11 clusters, and subsequent differential gene expression is obtained using t-test across Louvain clusters with Benjamini-Hochberg correction. We define differentially expressed genes (DEGs) using a 5% adjusted p-value cutoff and rank significant genes by the magnitude of their z-scores. We show in Figure 4a the data in PC-space colored with the respective methods, notably hierarchical clustering for Cluster LOCO-MP, neuralized KMeans for LRP and Louvain clustering for differential gene expression where labels have been aligned with the reported purified labels. Cluster LOCO-MP was also run with B = 5000 minipatches and minipatch ratio sizes αN = 0.42 for observations, αM = 0.26 for features for the cluster-level analysis. The cluster-level feature importance are here analyzed for three well separated cell types (CD14+ monocytes, dendritic cells, CD56+ natural killer cells) and one cell type that is harder to cluster (CD34+ hematopoietic stem cells) shown in the PC space in Figure 4a. 15

Reported purified labels

b.

Cluster LOCO-MP (HC)

Cluster LOCO-RAMPART (HC)

LYZ CST3 HLA-DRA HLA-DQA1 CD74 HLA-DRB1 AIF1 FTL HLA-DPB1 LST1

PC2

a.

PC1

Hierarchical Clustering

0.0

LYZ FTL AIF1 HLA-DRA LST1 CST3 CD74 HLA-DRB1 FCGR3A * HLA-DQA1 0.2

0.4

0.6

0.8

1.0

0.0

0.2

PC2

LRP (neuralized KMeans)

PC1

KMeans Clustering

PC2 CD4+/CD25 T Reg CD4+/CD45RA+/CD25- Naive T CD4+/CD45RO+ Memory CD8+ Cytotoxic T CD8+/CD45RA+ Naive Cytotoxic

CD14+ Monocyte CD19+ B CD34+ CD56+ NK Dendritic

0.2

0.4

0.6

0.8

1.0

1.0

0.0

0.2

0.4

0.6

0.8

1.0

0.8

1.0

PFI (KMeans)

RINT1 GCHFR NDUFB10 DAAM1 NUCB2 ISG20 STMN1 POU2AF1 BFAR TMEM256 0.0

0.8

LYZ IGLL1 GATA2 FTL CST3 EGFL7 RINT1 FAM185A AIF1 RP3-467N11.1

PBFI (KMeans)

PC1

0.6

IMPACC (HC)

FTL HLA-DRA LST1 CD74 AIF1 HLA-DQA1 LYZ PSAP HLA-DRB1 TYROBP * 0.0

0.4

HES4 ADRBK1 DPEP2 SLC11A1 KLF6 ETFA CFD GPX4 LRRC25 FGFR1OP2 0.2

0.4

0.6

0.8

1.0

0.0

0.2

0.4

0.6

Figure 3: Top 10 most globally important genes obtained via different feature importance methods. a. Normalized gene expression, reported and clustered via base clustering models - Hierarchical Clustering and standard KMeans clustering, in PC space. b. Feature importance scores: we color the genes by known marker genes extracted from the Azimuth Human PBMC database Hao et al. [19] and the manual annotations from Ding et al. [9] and Stuart et al. [53]. Genes that appear not to be in any known PBMC cell type marker set are colored in grey. We note that Cluster LOCO-MP finds the most marker genes, KMeans-based methods fail to identify known markers in the important genes.

16

Cluster LOCO: Feature Importance for Interpreting Clusters

Hierarchical Clustering

6

1 46 7

7

9

PC2

PC2

1 402 3 5

8

5

0.5

c.

CD14+ Monocyte (Cluster 5)

1.0

Dendritic (Cluster 9)

FTL LST1 AIF1 HLA-DRA PSAP FCGR3A * CTSS TYROBP * CD74 HLA-DQA1

1.0

0.0

LRP (KMeans)

IRF8 HLA-DQA1 AMICA1 VIM FCER1A CLEC10A LAT2 HLA-DQB1 HLA-DRA S100A10

0.5

1.0

0.0

1.0

Dendritic (Cluster 9)

(Cluster 5)

FCGR3A * FTL AIF1 LST1 PSAP FCER1G * TYROBP * CTSS CFD SPI1

0.5

1.0

0.0

CD56+ NK (Cluster 8)

1.0

0.0

1.0

0.0

1.0

SERPINB1 SNHG7 C19ORF77 IMPDH2 NPM1 HNRNPA1 HSP90AB1 RPS24 NFE2 C6ORF48

0.5

1.0

0.0

0.5

1.0

CD34+ (Cluster 7)

NKG7 * CTSW GZMA GNLY * PTPRCAP CD7 CST7 * GZMM GZMB CCL5 *

0.5

0.5

CD34+ (Cluster 7)

(Cluster 8)

LYZ CST3 HLA-DQA1 HLA-DQA2 HLA-DRA HLA-DPB1 HLA-DPA1 HLA-DRB1 FCER1A CD74

8: CD56+ NK 9: Dendritic

C19ORF77 PRSS57 NFE2 EGFL7 SERPINB1 HNRNPA1 EBPL CYTL1 EGR1 SNHG7

Louvain Clustering + DEG CD56+ NK

d. CD14+ Monocyte

0.5

0.0

5

CD34+ (Cluster 7)

NKG7 * GNLY * CTSW GZMB GZMA GZMH PRF1 CD7 CD63 IGFBP7

0.5

6: CD19+ B 7: CD34+

(Cluster 8)

NKG7 * SPON2 GZMB IGFBP7 GNLY * CLIC3 PRF1 CD7 S1PR5 VIM

0.5

3

PC1

Cluster LOCO-MP (HC)CD56+ NK

0.0

9 4

8

5

4: CD8+/CD45RA+ Naive Cytotoxic 5: CD14+ Monocyte

LYZ CD74 CST3 HLA-DRA FTL CD79A * FCGR3A * HLA-DRB1 CFD AIF1

10

6 7

2 0

PC1

Dendritic (Cluster 9)

(Cluster 5)

0.0

1

94 2

6 7

8

2: CD4+/CD45RO+ Memory 3: CD8+ Cytotoxic T

b. CD14+ Monocyte FTL LST1 AIF1 LYZ CTSS FCGR3A * HLA-DQA1 TYROBP * CFD FCER1G *

Louvain clustering

03

PC1

0: CD4+/CD25 T Reg 1: CD4+/CD45RA+/CD25- Naive T

0.0

1

9

30 8

PC1

0.0

2

KMeans clustering

PC2

Reported labels

PC2

a.

HNRNPA1 NPM1 SNHG7 C6ORF48 LDHB * EIF3L SNHG8 RPS24 IMPDH2 H2AFY

0.5

1.0

0.0

0.5

1.0

Figure 4: Top 10 most important genes for clusters corresponding to well separated cell types: CD14+ Monocytes (Cluster 5), Dendritic cells (Cluster 9), Natural Killer cells (Cluster 8) and CD34+ hematopoietic stem cells (Cluster 7). a. Normalized gene expression in PC space colored by cluster labels aligned with the reported labels. KMeans labels are from the neuralized KMeans model. b. Cluster LOCO-MP feature importance scores: Hierarchical Clustering was used as base clustering algorithm. c. Layer-wise Relevant Propagation (LRP) feature importance scores for neuralized KMeans. d. Louvain clustering and differential gene expression: genes ranked in size effect among those tested significant at level 0.05 in differential expression testing. Our method Cluster LOCOMP captures more consistent cluster-specific markers. 17

For each aligned cluster, we compare the top-ranked genes from Cluster LOCO-MP, LRP, and the standard Scanpy workflow (Figure 4b–d), with full top-10 rankings for all celltype-aligned clusters reported in the Additional figures (Figures 10, 11, 12). Across all clusters, Cluster LOCO-MP recovers a larger number of cluster-specific marker genes than the competing methods. In particular, for the well separated clusters, Cluster LOCOMP and LRP produce broadly similar rankings, and both are largely consistent with the DEG-based rankings. However, LRP and DEGs rank several genes with weaker cell-type specificity or genes absent from our marker reference set (CFD, VIM, LAT2, SPI1, S100A10, PTPRCAP ). We also note that Cluster LOCO-MP and LRP scores for monocytes and dendritic reflect the mixing of true dendritic cells with monocytes as shown in the reported labels of Figure 3a, with the presence of shared markers. On the other hand, this cell type is largely over-clustered via Louvain (see Figure 12) where the DEGs of cluster 4, 9 and the unmatched cluster 10 are very similar. In this particular case study, DEGs show possible false positive signal, with larger amount of markers significant for clusters that do not match the cell type they identify (see table 2). For the CD34+ hematopoietic stem cell cluster, which is less cleanly separated than the monocyte, dendritic-cell, and NK clusters, Cluster LOCO-MP is largely more interpretable than LRP and differential gene expression analysis. In this case, Cluster LOCO-MP ranks three known CD34+ HSC markers among its top 10 genes, namely PRSS57, EGFL7, and CYTL1, whereas LRP and differential gene expression rank no known CD34+ HSC markers, in fact LDHB appears as a DEG although reported to rather be a marker for Naive T cells or Memory T cells. This pattern is consistent with the geometry of the learned clustering shown in the PC space: well-separated clusters yield more specific and interpretable markers, while less cleanly separated populations, particularly Tcell subtypes, lead to more ambiguous marker rankings across all methods. We note however that this validation against known biological markers is necessarily incomplete: genes absent from a reference marker set may reflect noise or poor specificity, but they may also represent missing markers, or previously under-characterized context-specific genes. Nevertheless, Cluster LOCO-MP most consistently prioritizes cluster-specific markers, suggesting that its local feature importance scores better reflect the features supporting the learned clustering structure.

5 Discussion Clustering is often subject to debate regarding its rigor: is clustering an art or science [38]? In this work, we contribute to making clustering more scientifically rigorous by addressing a gap in the interpretability of clustering solutions. We devised Cluster LOCO, a novel family of feature importance metrics for clustering that are model-agnostic, flexible and scalable. Through synthetic and real-world evaluations, we have shown our method’s flexibility and Cluster LOCO provides feature-level explanation improving the reliability and trustworthiness of clustering solutions, especially for downstream analysis. In particular, in our single-cell transcriptomics application, we addressed a current challenge in single-cell analysis by providing a new unified framework for gene marker annotation that informs practitioners with interpretable, cluster-specific genes corresponding to cell-types when the clustering algorithm identifies a cluster aligned to the biological cell states.

18

Cluster LOCO: Feature Importance for Interpreting Clusters

Furthermore, Cluster LOCO’s model-agnostic design extends beyond typical clustering algorithms to encompass entire clustering pipelines: our framework can be applied on workflows that incorporate preliminary dimensionality reduction, or truly blackbox models, provided they can be applied to obtain new clustering solutions. Cluster LOCO was therefore implemented as an open-source Python library compatible with scikit-learn-style estimators and requires standard fit and predict functionality from the clustering model. As such, Cluster LOCO is universally applicable to both simple algorithms and sophisticated unsupervised architectures, while retaining a straightforward feature interpretation. However, several limitations remain: first, although the minipatch formulation improves scalability, Cluster LOCO still requires repeated clustering because of feature occlusion, making it more computationally demanding than simpler post-hoc scores. For instance, in settings involving neural network-based clustering models, architecture-specific interpretability methods may be more computationally efficient, although typically less general than our proposed approach. Second, Cluster LOCO focuses on feature-level interpretability: it explains which variables are important for a clustering solution globally or cluster-wise, but it does not directly provide sample-level local explanations for why a particular observation was assigned to a particular cluster. Future work could extend Cluster LOCO along both methodological and computational directions. Computationally, future work could further improve the efficiency of the minipatch procedure through adaptive sampling designs (e.g. multi-armed bandits, or active learning methods). Methodologically, the notion of generalizability could be adapted to unsupervised tasks beyond clustering, including dimensionality reduction, topic modeling, or representation learning. In each case, LOCO-style feature importance would ask which features are necessary for the learned structure to remain stable and generalizable under refitting, providing a path toward feature-level interpretability for a broader class of unsupervised workflows. Finally, an important direction is to further investigate the use of Cluster LOCO in single-cell transcriptomics. Our immune-cell analysis provides one case study in which cell-type labels were available from external validation experiments, but future work could evaluate the framework across additional datasets, cell-type annotation settings, pre-processing pipelines, and biological contexts. In particular, it would be useful to study how Cluster LOCO compares with marker discovery pipelines based on referenceatlas cell-type classification. More broadly, we believe Cluster LOCO can support a more transparent and statistically grounded use of clustering in scientific discovery.

19

Code availability Code for Cluster LOCO and scripts for reproducing the main results accompanied with data accession instructions are available at https://github.com/DataSlingers/ClusterLOCO. The PBMC 68k dataset is publicly available under Scanpy datasets with additional data of the original 10X Genomics analysis. Acknowledgements

Appendix 5.1 Simulation details 5.1.1 Interlaced moons toy example The example is constructed from 2 interlaced moons (circular arcs) in both features 1, 2 and 3. We parametrize for the 3 clusters respectively: • Cluster 1: X1 = 0.35 + 1.35 sin(T1 ) + ϵx and X2 = 0.5 + 1.35 p cos(T0 ) + ϵy where T1 ∼ 2 U (0.05π, 1.05π) where ϵx , ϵy ∼ N (0, 0.08 ). Then let R = (X2 − 0.1)2 + (X1 + 0.05)2 , X̃3 = −0.8 sin(1.5R) + 0.2 log(|1 + X2 |) + 0.9 + ϵz , ϵz ∼ N (0, 0.082 ) as well. • Cluster 2: X1 = −0.05 + sin(T2 ) + ϵx and X2 = −0.05p + cos(T1 ) + ϵy where T2 ∼ 2 U (0.95π, 1.95π) where ϵx , ϵy ∼ N (0, 0.08 ). Then let R = (X2 − 0.1)2 + (X1 + 0.05)2 , X̃3 = 0.8 sin(1.5R) + 0.2 log(|1 + X2 |) − 0.7 + ϵz , ϵz ∼ N (0, 0.082 ) as well. • Cluster 3: X1 = −0.05 + 0.7 sin(Tp 3 ) + ϵx and X2 = 0.55 + 0.7 cos(T3 ) + ϵy with T3 ∼ U (−0.15π, 1.15π). Then let R = (X2 − 0.1)2 + (X1 + 0.05)2 , X̃3 = 0.8 sin(1.5R) + 0.5 log(|1 + X2 |) + 0.66 + ϵz , ϵz ∼ N (0, 0.082 ) as well. By construction, X̃3 is a non linear transformation of X1 and X2 and is particularly correlated with X2 . We decorrelate X̃3 from X2 to get the final X3 feature: X3 = X̃3 − Cov(X2 ,X̃3 V ar(X2 ) X2 . The data is then augmented with 2 features X4 , X5 ∼ U ([−1, 2]) independently. 5.1.2 Comparative simulations For each simulation setting, we generated labeled data (Xi , Yi )N i=1 where Yi ∈ {1, . . . , K} denotes the cluster label and each cluster k contains nk observations. The data generating process first produces a low-dimensional signal representation in latent dimension X̃ ∈ Rd0 and is embedded into higher-dimensional signal feature space Rd . We then append pure noise features, so Xi ∈ RM where M = d + dnoise . Gaussian Mixture. The standard Gaussian mixture model outputs generates the signal latent X̃i ∈ Rd0 , Yi ∈ {1, . . . , K}. Let α ∈ R>0 be a parameter controlling cluster separation, µk ∈ Rd0 , Σk ∈ Rd0 ×d0 , then observations for cluster k are sampled from X̃i |Yi = k ∼ N (αµk , Σk ) where Σk is obtained via an onion covariance structure [47] that builds a correlation matrix layer-by-layer, hence its name, and µk ’s coordinates are sampled uniformly at random between [−1, 1]. 20

Cluster LOCO: Feature Importance for Interpreting Clusters

Gamma Mixture. The Gamma mixture introduces non-Gaussian marginal distributions while allowing dependence among features through a Gaussian copula. For each cluster k, let Rk denote a cluster-specific correlation matrix and let Fkj be the CDF of a gamma distribution with shape parameter aj and cluster-specific scale parameter skj , Γ(aj , skj ) . We define the joint distribution of X̃i | Yi = k by the Gaussian copula
CRk (u1 , . . . , ud0 ) = ΦRk Φ−1 (u1 ), . . . , Φ−1 (ud0 ) , where ΦRk is the CDF of a N (0, Rk ) with Rk ∈ Rd0 ,d0 , and Φ is the CDF of N (0, I). The samples are generated by drawing Zi ∼ N (0, Rk ), setting Uij = Φ(Zij ), and then applying −1 the inverse gamma CDF X̃ij = Fkj (Uij ). This construction gives each coordinate a gamma marginal distribution while using Rk to control the dependence structure within cluster k. Interlaced Moons/Circles. We generate half moons and circles in a two dimensional lower dimension first using mixtures of circles and half-moons. For each cluster k, a radius rk and center ck ∈ R2 are sampled, and the cluster shape is chosen to be either a full circle or a half-circle from a given probability vector. Points are sampled along the corresponding circular arc with additive Gaussian noise: Xi = ck + rk (cos θi , sin θi ) + εi , where θi is sampled either on [0, 2π] for circles or on an interval of length π for half-moons, with a random orientation. Cluster centers are placed so that the pairwise overlap between the underlying circles is bounded by a fixed maximum overlap percentage. This produces nonlinear, partially interlaced cluster structures. Embedding data in higher dimensions and noise features For high-dimensional simulations, we embed the low-dimensional signal into a larger ambient space using random feature maps: we generate random orthonormal projection maps, bringing the X̃ ∈ RN ×d0 latent observations into P̂ X̃ = X ∈ RN ×d . We also construct cluster-specific embeddings, allowing different clusters to be transformed through different random maps. Finally, we add pure noise signal S ∈ RN ×dnoise sampled from distributions such as Gaussian, Gamma with a = 1, s = 1, Student-t, uniform, triangular, or Laplace distributions. The final data generated by our simulator becomes X = [X̂, S] ∈ RN ×M .

21

5.2 Additional figures

Figure 5: Comparison of Cluster LOCO-Split scores with different measures of error/dissimilarity averaged over 20 runs with standard errors for the toy example described in Section 3.1. Mean squared error and Hamming distance show the more difference in normalized scores compared to other feature importance scores but we see that this is an effect of the scaling, feature 5 has a "worse" effect for these two errors relative to the rest of the scores, normalizing inflates the contribution of feature 2, 3 and 4 since all scores are normalized to be positive with highest importance at 1. Cluster LOCO

Stability LOCO

Silhouette LOCO

0.00 0.05 0.10 Permutation FI

0.2 0.0 Permutation FI

0.00 0.02 0.04 Fuzzy c-Means SHAP

Cluster LOCO

Stability LOCO

0.1 0.2 NEON (LRP)

0.0 0.2 Prototype FI

Silhouette LOCO

Benchmark methods

LOCO scores

X1 X2 X3 X4 X5 0.0

0.00

0.05 Method Spectral KMeans

X1 X2 X3

Feature type Signal Noise

X4 X5 0.0

0.5

0.0

0.2

0.4 0.000 0.005 0.010

0.00

0.05

0.10

0.0

0.2

ARIKMeans = 0.40 | ARISpectral = 1.00

Figure 6: Extended comparison of Cluster LOCO-Split with existing feature importance scores: we compared other LOCO-style scores with stability metric and silhouette score, for both KMeans and Spectral Clustering.

22

Cluster LOCO: Feature Importance for Interpreting Clusters

B. Gamma Mixture

C. Moons/Circles

1.0

1.0

1.0

0.8

0.8

0.8

0.6 0.4

0.6

Top-k Hits

Top-k Hits

Top-k Hits

A. Gaussian Mixture

0.4 0.2

0.2

0.6 0.4 0.2

0.0 20

50

200

Number of Features p

500

1000

20

200

Number of Features p

500

1000

1.0

1.0

0.4

0.4 0.2

0.2

50

20

50

200

500

1000

200

500

1000

Number of Features p

0.8

0.6

Top-k Hits

Top-k Hits

0.6

20

1.0

0.8

0.8

Top-k Hits

50

0.6 0.4 0.2

0.0 20

50

200

Number of Features p

Cluster LOCO RAMPART

500

1000

Cluster LOCO-MP

20

50

200

Number of Features p

PBFI (KMeans)

500

LRP (KMeans)

1000

IMPACC

Number of Features p

Permutation

Fuzzy C-Means SHAP

Figure 7: Top-k hits in three simulation settings for signal features p∗ = 10 and noise features features pnoise ∈ {10, 40, 90, 490, 990} with K = 3 clusters with N = 500 observations per cluster. Cluster LOCO-MP was run with αM = αN = 0.2 and B = 5000, Cluster LOCO-RAMPART with Brampart = 1000. Model agnostic methods are obtained with KMeans for Gaussian mixtures, Gamma mixture EM for Gamma mixture, Spectral Clustering for Moons and Circles. Model-specific methods are reported with their base model in the legend. Top-k hits were reported averaged over 100 replicates.

23

B. Gamma Mixture

C. Moons/Circles

1.0

1.0

1.0

0.8

0.8

0.8

0.6 0.4

0.6

Top-k Hits

Top-k Hits

Top-k Hits

A. Gaussian Mixture

0.4 0.2

0.2

0.6 0.4 0.2

0.0 20

50

200

Number of Features p

500

1000

20

50

200

Number of Features p

500

1000

50

20

50

200

500

1000

200

500

1000

Number of Features p

1.0

1.0 0.8

0.8

0.6 0.4

0.6

Top-k Hits

Top-k Hits

0.8

Top-k Hits

20

0.4 0.2 0.0

50

200

Number of Features p

Cluster LOCO RAMPART

500

1000

Cluster LOCO-MP

0.4 0.2

0.2 20

0.6

20

50

200

Number of Features p

PBFI (KMeans)

500

LRP (KMeans)

1000

IMPACC

0.0

Number of Features p

Permutation

Fuzzy C-Means SHAP

Figure 8: Top-k hits in three simulation settings for signal features p∗ = 10 and noise features features pnoise ∈ {10, 40, 90, 490, 990} with K = 3 clusters with N = 300 observations per cluster. Cluster LOCO-MP was run with B = 5000, Cluster LOCO-RAMPART with Brampart = 1000 and adaptive αN ∈ (0.1, 0.5), αp ∈ (0.1, 0.5). Model agnostic methods are obtained with KMeans for Gaussian mixtures, Gamma mixture EM for Gamma mixture, Spectral Clustering for Moons and Circles. Model-specific methods are reported with their base model in the legend. Top-k hits were reported averaged over 100 replicates.

24

Cluster LOCO: Feature Importance for Interpreting Clusters

B. Gamma Mixture

C. Moons/Circles

1.0

1.0

0.8

0.8

0.8

0.6

0.6

0.6

0.4

0.4 0.2

0.2 0.0

Top-k Hits

1.0

Top-k Hits

Top-k Hits

A. Gaussian Mixture

0.2

0.0 20

50

200

Number of Features p

500

1000

0.4

0.0 20

50

200

Number of Features p

500

1000

20

50

20

50

200

500

1000

200

500

1000

Number of Features p

1.0

1.0 0.8

0.8

0.4

0.4 0.2

0.2 0.0

0.6

Top-k Hits

0.6

Top-k Hits

Top-k Hits

0.8

50

200

Number of Features p

Cluster LOCO RAMPART

500

1000

Cluster LOCO-MP

0.4 0.2

0.0 20

0.6

0.0 20

50

200

Number of Features p

PBFI (KMeans)

500

LRP (KMeans)

1000

IMPACC

Number of Features p

Permutation

Fuzzy C-Means SHAP

Figure 9: Top-k hits in three simulation settings for signal features p∗ = 10 and noise features features pnoise ∈ {10, 40, 90, 490, 990} with K = 3 clusters with N = 100 observations per cluster. Cluster LOCO-MP was run with αM = αN = 0.2 and B = 5000, Cluster LOCO-RAMPART with Brampart = 1000. Model agnostic methods are obtained with KMeans for Gaussian mixtures, Gamma mixture EM for Gamma mixture, Spectral Clustering for Moons and Circles. Model-specific methods are reported with their base model in the legend. Top-k hits were reported averaged over 100 replicates.

25

CD14+ Monocyte (Cluster 5) FTL LST1 AIF1 LYZ CTSS FCGR3A * HLA-DQA1 TYROBP * CFD FCER1G *

0.0

0.5

1.0

NKG7 * SPON2 GZMB IGFBP7 GNLY * CLIC3 PRF1 CD7 S1PR5 VIM

0.0

CD34+ (Cluster 7)

0.5

1.0

1.0

0.0

0.5

1.0

0.5

PPIB MZB1 TNFRSF17 FKBP11 SUB1 UBE2J1 ISG20 AL928768.3 RABAC1 SPCS2

0.0

1.0

0.5

1.0

CD4+/CD25 T Reg (Cluster 0)

CD19+ B (Cluster 6) CD79B * CD79A * MS4A1 * CD37 BLK LTB SMARCB1 LYZ BANK1 CD74

0.0

CD4+/CD45RA+/CD25- Naive T (Cluster 1)

GNLY * CCL5 * NKG7 * GZMA CTSW GZMH HOPX PRF1 CST7 * CD7

Dendritic (Cluster 9) LYZ CD74 CST3 HLA-DRA FTL CD79A * FCGR3A * HLA-DRB1 CFD AIF1

0.5

CD8+ Cytotoxic T (Cluster 3)

CD56+ NK (Cluster 8)

RP11-291B21.2 ITM2A CNN2 CD8B * C9ORF78 PQBP1 TPM4 SKP1 CORO1B CD52

C19ORF77 PRSS57 NFE2 EGFL7 SERPINB1 HNRNPA1 EBPL CYTL1 EGR1 SNHG7

0.0

Cluster LOCO-MP (HC)

CD4+/CD45RO+ Memory (Cluster 2)

0.5

1.0

0.5

1.0

CD8+/CD45RA+ Naive Cytotoxic (Cluster 4)

HMGB2 STMN1 HLA-DQA1 HLA-DPA1 HLA-DRB1 LAMTOR1 HLA-DRB5 F12 EIF4A1 ST6GALNAC1

0.0

0.0

IRF8 CPVL CLEC10A FCER1A PBXIP1 NME4 IL32 UFD1L MIEN1 RNF130

0.0

0.5

1.0

0.0

0.5

1.0

Figure 10: Top 10 most important genes for clusters via Cluster LOCO-MP feature importance scores: Hierarchical Clustering was used as base clustering algorithm. Labels were aligned with the reported purified labels. Known markers from our reference set are colored with the cell-type they identify.

CD4+/CD45RA+/CD25- Naive T (Cluster 1) PPIB MZB1 SPCS2 TNFRSF17 TMBIM6 UBE2J1 MANF SSR4 FKBP11 LTB

0.0

DEK LAT2 CD33 ATP5I DHRS7 RNH1 GSTP1 DDX5 H2AFY EIF4A1

0.5

1.0

0.0

CD4+/CD25 T Reg (Cluster 0)

0.5

1.0

0.0

1.0

0.0

0.5

1.0

0.0

1.0

0.0

FTL LST1 AIF1 HLA-DRA PSAP FCGR3A * CTSS TYROBP * CD74 HLA-DQA1

0.5

1.0

1.0

0.0

0.5

1.0

CD4+/CD45RO+ Memory (Cluster 2)

JUNB TPM4 LTB TRAPPC1 DDX5 PPIA CALM3 BTG1 STMN1 TXNIP

0.5

0.0

CD8+ Cytotoxic T (Cluster 3)

SERPINB1 SNHG7 C19ORF77 IMPDH2 NPM1 HNRNPA1 HSP90AB1 RPS24 NFE2 C6ORF48

0.5

CD14+ Monocyte (Cluster 5)

Dendritic (Cluster 9) IRF8 HLA-DQA1 AMICA1 VIM FCER1A CLEC10A LAT2 HLA-DQB1 HLA-DRA S100A10

CD34+ (Cluster 7)

CD79B * CD74 MS4A1 * LTB CST3 CD79A * LYZ HLA-DRA CD37 PTPRCAP

0.5

CD56+ NK (Cluster 8)

NKG7 * GNLY * CTSW GZMB GZMA GZMH PRF1 CD7 CD63 IGFBP7

CD19+ B (Cluster 6)

TPM4 TRAPPC1 CNN2 CALM3 TXNIP JUNB DDX5 LTB ATP5B PPIA

0.0

LRP (KMeans)

CD8+/CD45RA+ Naive Cytotoxic (Cluster 4)

CD33 HLA-DQA1 DEK FCER1A CLEC10A HLA-DQA2 LGALS2 HLA-DRA RNF130 HLA-DQB1

0.5

1.0

0.0

0.5

1.0

Figure 11: Top 10 most important genes for clusters via LRP feature importance scores, LRP uses a neuralized KMeans clustering algorithm. Labels were aligned with the reported purified labels. Known markers from our reference set are colored with the cell-type they identify.

26

Cluster LOCO: Feature Importance for Interpreting Clusters

Louvain Clustering + DEG

CD4+/CD25 T Reg (cluster 0)

CD14+ Monocyte (cluster 5)

CD3D * LDHB * AES CD3E * GIMAP7 PTPRCAP CD52 IL32 NOSIP * LTB

Dendritic (cluster 9)

FCGR3A * FTL AIF1 LST1 PSAP FCER1G * TYROBP * CTSS CFD SPI1

0.0

0.2

0.4

0.6

0.8

1.0

LYZ CST3 HLA-DQA1 HLA-DQA2 HLA-DRA HLA-DPB1 HLA-DPA1 HLA-DRB1 FCER1A CD74

0.0

0.2

0.4

0.6

0.8

1.0

CD8+ Cytotoxic T (cluster 3)

CD19+ B (cluster 6) CD79A * CD79B * PTPRCAP MS4A1 * CD74 CD52 HLA-DPB1 HLA-DRA LTB CD37 0.2

0.4

0.6

0.8

1.0

0.2

0.4

0.6

0.8

1.0

0.2

0.4

0.6

0.8

1.0

0.4

0.6

0.8

1.0

0.0

0.0

0.0

0.2

0.4

0.4

0.6

0.8

1.0

0.2

0.4

0.6

0.8

1.0

0.8

1.0

0.0

0.2

0.4

CD34+ (cluster 7)

0.0

0.2

0.4

0.6

Figure 12: Top 10 most important DEGs obtained with Louvain clustering and t-test with BH correction. DEGs were filtered at significance level of 0.05 and ranked by z-score magnitude. Labels were aligned with the reported purified labels. Known markers from our reference set are colored with the cell-type they identify.

27

0.8

1.0

Unmatched cluster

HNRNPA1 NPM1 SNHG7 C6ORF48 LDHB * EIF3L SNHG8 RPS24 IMPDH2 H2AFY 0.2

0.6

CST3 LYZ HLA-DPA1 HLA-DQA1 HLA-DRA HLA-DPB1 CD74 HLA-DRB1 CPVL IRF8

CD4+/CD45RO+ Memory (cluster 2) STMN1 HMGB2 CALM3 CD52 PPIA HNRNPA1 CD3D * PTPRCAP CALM1 LDHB *

0.0

0.2

LYZ CST3 HLA-DRA HLA-DRB1 CD74 HLA-DPB1 FCER1A HLA-DRB5 HLA-DQA1 S100A9

0.0

CD4+/CD45RA+/CD25- Naive T (cluster 1) PPIB MZB1 FKBP11 SPCS2 ISG20 IGJ SUB1 TNFRSF17 TMBIM6 SSR4

0.0

NKG7 * CTSW GZMA GNLY * PTPRCAP CD7 CST7 * GZMM GZMB CCL5 *

CD8+/CD45RA+ Naive Cytotoxic (cluster 4)

FTL LST1 AIF1 TYROBP * FCER1G * FCGR3A * CTSS CFD PSAP PILRA

0.0

CD56+ NK (cluster 8)

0.6

0.8

1.0

5.3 Additional tables Algorithm 3 Cluster LOCO-RAMPART Require: Unlabeled data X ∈ RN ×M , top-k target k, minipatch sizes n, m, minipatches per round B. 1: Set t = 0. Initialize active feature set S0 = [M ]. 2: while |St | > k do 3: Obtain round-t Cluster LOCO-MP estimates: (t)

ˆ }j∈St ← Cluster LOCO-MP(X·,St , n, m, B) {∆ j t ) Determine Cluster LOCO-MP estimates’ ranks (in ascending order): (r̂1t , . . . , r|S t| t t at feature (τ̂1 , . . . , τ̂|St | ) respectively t 5: Retain the top half of the candidate set St+1 → {τ̂1t , · · · , τ̂|S }. t |/2 6: Set t ← t + 1. 7: end while ˆ (T ) }j∈S , |ST | ≈ k. Ensure: {∆ T j

4:

Cluster LOCO-MP LRP DEGs

Cluster 0

Cluster 1

Cluster 2

Cluster 3

Cluster 4

Cluster 5

Cluster 6

Cluster 7

Cluster 8

Cluster 9

CD4+ T reg

CD4+ Naive T

CD4+ Mem T

CD8+ Cytotoxic T

CD8+ Naive T

CD14+ Monocyte

CD19+ B

CD34+

CD56+ NK

Dendritic

0/10

0/10

0/10

5/10

0/10

5/10

4/10

3/10

5/10

5/10

0/10 0/10

0/10 0/10

0/10 1/10

0/10 1/10

0/10 0/10

4/10 4/10

3/10 3/10

0/10 0/10

4/10 4/10

5/10 8/10

Table 1: Number of markers in the top-10 genes matching known markers of identified cell-type. We report in blue the highest proportion of consistent marker identified across methods for each cell types.

Cluster LOCO-MP LRP DEGs

Cluster 0

Cluster 1

Cluster 2

Cluster 3

Cluster 4

Cluster 5

Cluster 6

Cluster 7

Cluster 8

Cluster 9

CD4+ T reg

CD4+ Naive T

CD4+ Mem T

CD8+ Cytotoxic T

CD8+ Naive T

CD14+ Monocyte

CD19+ B

CD34+

CD56+ NK

Dendritic

3/10

1/10

1/10

3/10

3/10

2/10

3/10

2/10

1/10

3/10

1/10 6/10

2/10 1/10

5/10 2/10

2/10 4/10

0/10 9/10

3/10 1/10

5/10 4/10

0/10 1/10

3/10 3/10

0/10 1/10

Table 2: Number of markers in the top-10 genes corresponding to known markers for other cell-types (possible false positive signal). We report in red the highest proportion of inconsistent markers across methods identified for each cell type.

28

Cluster LOCO: Feature Importance for Interpreting Clusters

References [1] Salem Alelyani, Jiliang Tang, and Huan Liu. Feature Selection for Clustering: A Review. In Data Clustering. Chapman and Hall/CRC, 2014. ISBN 978-1-315-37351-5. Num Pages: 32. [2] Genevera I. Allen, Luqin Gan, and Lili Zheng. Interpretable Machine Learning for Discovery: Statistical Challenges \& Opportunities, August 2023. URL http://arxiv. org/abs/2308.01475. arXiv:2308.01475 [stat.ML]. [3] Asa Ben-Hur, Andre Elisseeff, and Isabelle Guyon. A stability based method for discovering structure in clustered data. In Biocomputing 2002, Kauai, Hawaii, USA, December 2001. WORLD SCIENTIFIC. doi: 10.1142/9789812799623_0002. URL http://www.worldscientific.com/doi/abs/10.1142/9789812799623_0002. [4] Leo Breiman. Random Forests. Machine Learning, 45(1):5–32, October 2001. ISSN 1573-0565. doi: 10.1023/A:1010933404324. URL https://doi.org/10.1023/A: 1010933404324. [5] Yuxi Chen, Tiffany Tang, and Genevera Allen. Top-$k$ Feature Importance Ranking, September 2025. URL http://arxiv.org/abs/2509.15420. arXiv:2509.15420 [cs]. [6] M. Dash, K. Choi, P. Scheuermann, and Huan Liu. Feature selection for clustering - a filter solution. In 2002 IEEE International Conference on Data Mining, 2002. Proceedings., pages 115–122, December 2002. doi: 10.1109/ICDM.2002.1183893. URL https://ieeexplore.ieee.org/abstract/document/1183893. [7] Alan DenAdel, Michelle L. Ramseier, Andrew W. Navia, Alex K. Shalek, Srivatsan Raghavan, Peter S. Winter, Ava P. Amini, and Lorin Crawford. A knockoff calibration method to avoid over-clustering in single-cell RNA-sequencing, March 2024. URL https://www.biorxiv.org/content/10.1101/2024.03.08.584180v1. Pages: 2024.03.08.584180 Section: New Results. [8] Chris Ding and Xiaofeng He. K -means clustering via principal component analysis. In Twenty-first international conference on Machine learning - ICML ’04, page 29, Banff, Alberta, Canada, 2004. ACM Press. doi: 10.1145/1015330.1015408. URL http: //portal.acm.org/citation.cfm?doid=1015330.1015408. [9] Jiarui Ding, Xian Adiconis, Sean K. Simmons, Monika S. Kowalczyk, Cynthia C. Hession, Nemanja D. Marjanovic, Travis K. Hughes, Marc H. Wadsworth, Tyler Burks, Lan T. Nguyen, John Y. H. Kwon, Boaz Barak, William Ge, Amanda J. Kedaigle, Shaina Carroll, Shuqiang Li, Nir Hacohen, Orit Rozenblatt-Rosen, Alex K. Shalek, Alexandra-Chloé Villani, Aviv Regev, and Joshua Z. Levin. Systematic comparison of single-cell and single-nucleus RNA-sequencing methods. Nature Biotechnology, 38 (6):737–746, June 2020. ISSN 1546-1696. doi: 10.1038/s41587-020-0465-8. URL https://www.nature.com/articles/s41587-020-0465-8. [10] Martin Ester, Hans-Peter Kriegel, and Xiaowei Xu. A Density-Based Algorithm for Discovering Clusters in Large Spatial Databases with Noise. AAAI, 1996. 29

[11] Pedro F. Ferreira, Jack Kuipers, and Niko Beerenwinkel. Identifying hierarchical cell states and gene signatures with deep exponential families for single-cell transcriptomics, January 2024. URL https://www.biorxiv.org/content/10.1101/2022.10. 15.512383v2. Pages: 2022.10.15.512383 Section: New Results. [12] Aaron Fisher, Cynthia Rudin, and Francesca Dominici. All Models are Wrong, but Many are Useful: Learning a Variable’s Importance by Studying an Entire Class of Prediction Models Simultaneously. Journal of machine learning research : JMLR, 20:177, 2019. ISSN 1532-4435. URL https://pmc.ncbi.nlm.nih.gov/articles/ PMC8323609/. [13] Luqin Gan and Genevera I. Allen. Fast and interpretable consensus clustering via minipatch learning. PLOS Computational Biology, 18(10):e1010577, October 2022. ISSN 1553-7358. doi: 10.1371/journal.pcbi.1010577. URL https://journals.plos.org/ ploscompbiol/article?id=10.1371/journal.pcbi.1010577. [14] Luqin Gan, Lili Zheng, and Genevera I. Allen. Model-Agnostic Confidence Intervals for Feature Importance: A Fast and Powerful Approach Using Minipatch Ensembles, January 2023. URL http://arxiv.org/abs/2206.02088. arXiv:2206.02088 [stat]. [15] Luqin Gan, Tarek M. Zikry, and Genevera I. Allen. Are machine learning interpretations reliable? A stability study on global interpretations, May 2025. URL http://arxiv. org/abs/2505.15728. arXiv:2505.15728 [stat]. [16] Greg Gibson. Perspectives on rigor and reproducibility in single cell genomics. PLoS Genetics, 18(5):e1010210, May 2022. ISSN 1553-7390. doi: 10.1371/journal.pgen.1010210. URL https://pmc.ncbi.nlm.nih.gov/articles/PMC9122178/. [17] Mark S. Handcock, Adrian E. Raftery, and Jeremy M. Tantrum. Model-Based Clustering for Social Networks. Journal of the Royal Statistical Society Series A: Statistics in Society, 170(2):301–354, March 2007. ISSN 0964-1998, 1467-985X. doi: 10.1111/j. 1467-985X.2007.00471.x. URL https://academic.oup.com/jrsssa/article/170/2/ 301/7085294. [18] Julia Handl, Joshua Knowles, and Douglas B. Kell. Computational cluster validation in post-genomic data analysis. Bioinformatics, 21(15):3201–3212, August 2005. ISSN 1367-4803. doi: 10.1093/bioinformatics/bti517. URL https://doi.org/10. 1093/bioinformatics/bti517. [19] Yuhan Hao, Stephanie Hao, Erica Andersen-Nissen, William M. Mauck, Shiwei Zheng, Andrew Butler, Maddie J. Lee, Aaron J. Wilk, Charlotte Darby, Michael Zager, Paul Hoffman, Marlon Stoeckius, Efthymia Papalexi, Eleni P. Mimitou, Jaison Jain, Avi Srivastava, Tim Stuart, Lamar M. Fleming, Bertrand Yeung, Angela J. Rogers, Juliana M. McElrath, Catherine A. Blish, Raphael Gottardo, Peter Smibert, and Rahul Satija. Integrated analysis of multimodal single-cell data. Cell, 184(13):3573– 3587.e29, June 2021. ISSN 0092-8674, 1097-4172. doi: 10.1016/j.cell.2021.04.048. URL https://www.cell.com/cell/abstract/S0092-8674(21)00583-3. 30

Cluster LOCO: Feature Importance for Interpreting Clusters

[20] Lianyu Hu, Mudi Jiang, Junjie Dong, Xinying Liu, and Zengyou He. Interpretable Clustering: A Survey, September 2024. URL http://arxiv.org/abs/2409.00743. arXiv:2409.00743 [cs]. [21] Jiashun Jin and Wanjie Wang. Influential features PCA for high dimensional clustering. The Annals of Statistics, 44(6):2323–2359, December 2016. ISSN 0090-5364, 2168-8966. doi: 10.1214/15-AOS1423. URL https://projecteuclid.org/journals/annals-of-statistics/volume-44/ issue-6/Influential-features-PCA-for-high-dimensional-clustering/10. 1214/15-AOS1423.full. [22] Ian T. Jolliffe and Jorge Cadima. Principal component analysis: a review and recent developments. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 374(2065):20150202, April 2016. ISSN 1364-503X. doi: 10.1098/rsta.2015.0202. URL https://doi.org/10.1098/rsta.2015.0202. [23] A. Karanikola, C. M. Liapis, and S. Kotsiantis. Investigating cluster validation metrics for optimal number of clusters determination. Intelligent Decision Technologies, 15(4): 809–824, 2021. [24] Jacob Kauffmann, Malte Esders, Lukas Ruff, Grégoire Montavon, Wojciech Samek, and Klaus-Robert Müller. From Clustering to Cluster Explanations via Neural Networks. IEEE Transactions on Neural Networks and Learning Systems, 35(2):1926–1940, February 2024. ISSN 2162-237X, 2162-2388. doi: 10.1109/TNNLS.2022.3185901. URL http://arxiv.org/abs/1906.07633. arXiv:1906.07633 [cs]. [25] Vladimir Yu Kiselev, Tallulah S. Andrews, and Martin Hemberg. Challenges in unsupervised clustering of single-cell RNA-seq data. Nature Reviews Genetics, 20(5):273– 282, May 2019. ISSN 1471-0056, 1471-0064. doi: 10.1038/s41576-018-0088-9. URL https://www.nature.com/articles/s41576-018-0088-9. [26] Hans-Peter Kriegel, Peer Kröger, and Arthur Zimek. Clustering high-dimensional data: A survey on subspace clustering, pattern-based clustering, and correlation clustering. ACM Trans. Knowl. Discov. Data, 3(1):1:1–1:58, March 2009. ISSN 1556-4681. doi: 10. 1145/1497577.1497578. URL https://dl.acm.org/doi/10.1145/1497577.1497578. [27] H. W. Kuhn. The Hungarian method for the assignment problem. Naval Research Logistics Quarterly, 2(1-2):83–97, 1955. ISSN 1931-9193. doi: 10.1002/nav.3800020109. URL https://onlinelibrary.wiley.com/doi/abs/10.1002/nav.3800020109. _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/nav.3800020109. [28] Tilman Lange, Volker Roth, Mikio L. Braun, and Joachim M. Buhmann. Stabilitybased validation of clustering solutions. Neural Computation, 16(6):1299–1323, June 2004. ISSN 0899-7667. doi: 10.1162/089976604773717621. [29] Jing Lei, Max G’Sell, Alessandro Rinaldo, Ryan J. Tibshirani, and Larry Wasserman. Distribution-Free Predictive Inference For Regression, March 2017. URL http://arxiv.org/abs/1604.04173. arXiv:1604.04173 [stat]. 31

[30] Diego De Leon-Oliva, Cielo Garcia-Montero, Oscar Fraile-Martinez, Diego Liviu Boaru, Luis García-Puente, Antonio Rios-Parra, Maria J. Garrido-Gil, Carlos CasanovaMartín, Natalio García-Honduvilla, Julia Bujan, Luis G. Guijarro, Melchor AlvarezMon, and Miguel A. Ortega. AIF1: Function and Connection with Inflammatory Diseases. Biology, 12(5):694, May 2023. doi: 10.3390/biology12050694. URL https://pmc.ncbi.nlm.nih.gov/articles/PMC10215110/. [31] E. Levine and E. Domany. Resampling method for unsupervised estimation of cluster validity. Neural Computation, 13(11):2573–2593, 2001. [32] Shenghao Li, Hui Guo, Simai Zhang, Yizhou Li, and Menglong Li. Attention-based deep clustering method for scRNA-seq cell type identification. PLOS Computational Biology, 19(11):e1011641, November 2023. ISSN 1553-7358. doi: 10.1371/journal.pcbi. 1011641. URL https://journals.plos.org/ploscompbiol/article?id=10.1371/ journal.pcbi.1011641. [33] Camille Little, Lili Zheng, and Genevera Allen. iLOCO: Distribution-Free Inference for Feature Interactions, May 2025. URL http://arxiv.org/abs/2502.06661. arXiv:2502.06661 [stat]. [34] Gilles Louppe. Understanding Random Forests: From Theory to Practice, June 2015. URL http://arxiv.org/abs/1407.7502. arXiv:1407.7502 [stat]. [35] Malte D. Luecken and Fabian J. Theis. Current best practices in single-cell RNAseq analysis: a tutorial. Molecular Systems Biology, 15(6):MSB188746, June 2019. ISSN 1744-4292. doi: 10.15252/msb.20188746. URL https://doi.org/10.15252/ msb.20188746. [36] Scott Lundberg and Su-In Lee. A Unified Approach to Interpreting Model Predictions, November 2017. URL http://arxiv.org/abs/1705.07874. arXiv:1705.07874 [cs]. [37] Ulrike von Luxburg. Clustering Stability: An Overview. Foundations and Trends® in Machine Learning, 2(3):235–274, 2009. ISSN 1935-8237, 1935-8245. doi: 10.1561/ 2200000008. URL http://arxiv.org/abs/1007.1075. arXiv:1007.1075 [stat]. [38] Ulrike von Luxburg, Robert C. Williamson, and Isabelle Guyon. Clustering: Science or Art? In Proceedings of ICML Workshop on Unsupervised and Transfer Learning, pages 65–79. JMLR Workshop and Conference Proceedings, June 2012. URL https: //proceedings.mlr.press/v27/luxburg12a.html. [39] David Lähnemann, Johannes Köster, Ewa Szczurek, Davis J. McCarthy, Stephanie C. Hicks, Mark D. Robinson, Catalina A. Vallejos, Kieran R. Campbell, Niko Beerenwinkel, Ahmed Mahfouz, Luca Pinello, Pavel Skums, Alexandros Stamatakis, Camille Stephan-Otto Attolini, Samuel Aparicio, Jasmijn Baaijens, Marleen Balvert, Buys de Barbanson, Antonio Cappuccio, Giacomo Corleone, Bas E. Dutilh, Maria Florescu, Victor Guryev, Rens Holmer, Katharina Jahn, Thamar Jessurun Lobo, Emma M. Keizer, Indu Khatri, Szymon M. Kielbasa, Jan O. Korbel, Alexey M. Kozlov, Tzu-Hao Kuo, Boudewijn P.F. Lelieveldt, Ion I. Mandoiu, John C. Marioni, Tobias Marschall, Felix Mölder, Amir Niknejad, Alicja Rączkowska, Marcel Reinders, Jeroen de Ridder, 32

Cluster LOCO: Feature Importance for Interpreting Clusters

Antoine-Emmanuel Saliba, Antonios Somarakis, Oliver Stegle, Fabian J. Theis, Huan Yang, Alex Zelikovsky, Alice C. McHardy, Benjamin J. Raphael, Sohrab P. Shah, and Alexander Schönhuth. Eleven grand challenges in single-cell data science. Genome Biology, 21(1):31, February 2020. ISSN 1474-760X. doi: 10.1186/s13059-020-1926-6. URL https://doi.org/10.1186/s13059-020-1926-6. [40] Masayoshi Mase, Art B. Owen, and Benjamin B. Seiler. Cohort Shapley value for algorithmic fairness. Technical Report arXiv:2105.07168, arXiv, May 2021. URL http: //arxiv.org/abs/2105.07168. arXiv:2105.07168 [cs, econ, stat] type: article. [41] J. Materne. The structure of nearby clusters of galaxies. Hierarchical clustering and an application to the Leo region. Astronomy and Astrophysics, 63:401–409, February 1978. ISSN 0004-6361. URL https://ui.adsabs.harvard.edu/abs/1978A&A....63..401M. ADS Bibcode: 1978A&A....63..401M. [42] Erxue Min, Xifeng Guo, Qiang Liu, Gen Zhang, Jianjing Cui, and Jun Long. A Survey of Clustering With Deep Learning: From the Perspective of Network Architecture. IEEE Access, 6:39501–39514, 2018. ISSN 2169-3536. doi: 10.1109/ACCESS.2018.2855437. URL https://ieeexplore.ieee.org/abstract/document/8412085. [43] Christoph Molnar. Interpretable Machine Learning (Third Edition). Leanpub, February 2018. URL https://leanpub.next/interpretable-machine-learning. [44] Gonzalo Nápoles, Niels Griffioen, Samaneh Khoshrou, and Çiçek Güven. Feature Importance for Clustering. In Verónica Vasconcelos, Inês Domingues, and Simão Paredes, editors, Progress in Pattern Recognition, Image Analysis, Computer Vision, and Applications, pages 31–45, Cham, 2024. Springer Nature Switzerland. ISBN 978-3-03149018-7. doi: 10.1007/978-3-031-49018-7_3. [45] Roy Oelen, Dylan H. de Vries, Harm Brugge, M. Grace Gordon, Martijn Vochteloo, Chun J. Ye, Harm-Jan Westra, Lude Franke, and Monique G. P. van der Wijst. Singlecell RNA-sequencing of peripheral blood mononuclear cells reveals widespread, contextspecific gene expression regulation upon pathogenic exposure. Nature Communications, 13(1):3267, June 2022. ISSN 2041-1723. doi: 10.1038/s41467-022-30893-5. URL https: //www.nature.com/articles/s41467-022-30893-5. [46] Oliver Pfaffel. FeatureImpCluster: Feature Importance for Partitional Clustering, May 2020. URL https://CRAN.R-project.org/package=FeatureImpCluster. Institution: Comprehensive R Archive Network Pages: 0.1.5. [47] Weiliang Qiu and Harry Joe. Generation of Random Clusters with Specified Degree of Separation. Journal of Classification, 23(2):315–334, September 2006. ISSN 1432-1343. doi: 10.1007/s00357-006-0018-y. URL https://doi.org/10.1007/ s00357-006-0018-y. [48] Alessandro Rinaldo, Ryan J. Tibshirani, and Larry Wasserman. Comment: Statistical Inference from a Predictive Perspective. Statistical Science, 34(4):599–603, 2019. doi: 10.1214/19-STS748. URL https://www.researchgate.net/publication/ 338469585_Comment_Statistical_Inference_from_a_Predictive_Perspective. 33

[49] Volker Roth and Tilman Lange. Feature Selection in Clustering Problems. In Advances in Neural Information Processing Systems, volume 16. MIT Press, 2003. URL https://proceedings.neurips.cc/paper_files/paper/2003/hash/ bb03e43ffe34eeb242a2ee4a4f125e56-Abstract.html. [50] Peter J. Rousseeuw. Silhouettes: A graphical aid to the interpretation and validation of cluster analysis. Journal of Computational and Applied Mathematics, 20: 53–65, November 1987. ISSN 0377-0427. doi: 10.1016/0377-0427(87)90125-7. URL https://www.sciencedirect.com/science/article/pii/0377042787901257. [51] Lloyd S. Shapley. A value for n-person games. In Contributions to the Theory of Games, volume 2, pages 307–317. Princeton University Press, 1953. [52] Dongyuan Song, Kexin Li, Xinzhou Ge, and Jingyi Jessica Li. ClusterDE: a postclustering differential expression (DE) method robust to false-positive inflation caused by double dipping. Research Square, page rs.3.rs, August 2023. doi: 10.21203/rs.3. rs-3211191/v1. URL https://pmc.ncbi.nlm.nih.gov/articles/PMC10418557/. [53] Tim Stuart, Andrew Butler, Paul Hoffman, Christoph Hafemeister, Efthymia Papalexi, William M. Mauck, Yuhan Hao, Marlon Stoeckius, Peter Smibert, and Rahul Satija. Comprehensive Integration of Single-Cell Data. Cell, 177(7):1888–1902.e21, June 2019. ISSN 0092-8674, 1097-4172. doi: 10.1016/j.cell.2019.05.031. URL https://www.cell. com/cell/abstract/S0092-8674(19)30559-8. [54] Mukund Sundararajan and Amir Najmi. The many Shapley values for model explanation. In Proceedings of the 37th International Conference on Machine Learning, volume 119 of ICML’20, pages 9269–9278. JMLR.org, July 2020. [55] Peter J. Thul and Cecilia Lindskog. The human protein atlas: A spatial map of the human proteome. Protein Science : A Publication of the Protein Society, 27(1):233– 244, January 2018. ISSN 0961-8368. doi: 10.1002/pro.3307. URL https://pmc.ncbi. nlm.nih.gov/articles/PMC5734309/. [56] Robert Tibshirani and Guenther Walther. Cluster Validation by Prediction Strength. Journal of Computational and Graphical Statistics, 14(3):511–528, September 2005. ISSN 1061-8600, 1537-2715. doi: 10.1198/106186005X59243. URL https://www. tandfonline.com/doi/full/10.1198/106186005X59243. [57] Mohammad Taha Toghani and Genevera I. Allen. MP-Boost: Minipatch Boosting via Adaptive Feature and Observation Sampling. In 2021 IEEE International Conference on Big Data and Smart Computing (BigComp), pages 75–78, January 2021. doi: 10.1109/BigComp51126.2021.00023. URL http://arxiv.org/abs/2011.07218. arXiv:2011.07218 [stat]. [58] Theresa Ullmann, Christian Hennig, and Anne-Laure Boulesteix. Validation of cluster analysis results on validation data: A systematic framework. WIREs Data Mining and Knowledge Discovery, 12(3):e1444, 2022. ISSN 1942-4795. doi: 10.1002/ widm.1444. URL https://onlinelibrary.wiley.com/doi/abs/10.1002/widm.1444. _eprint: https://wires.onlinelibrary.wiley.com/doi/pdf/10.1002/widm.1444. 34

Cluster LOCO: Feature Importance for Interpreting Clusters

[59] Isabella Verdinelli and Larry Wasserman. Feature Importance: A Closer Look at Shapley Values and LOCO, March 2023. URL http://arxiv.org/abs/2303.05981. arXiv:2303.05981 [stat]. [60] Alexandra-Chloé Villani, Rahul Satija, Gary Reynolds, Siranush Sarkizova, Karthik Shekhar, James Fletcher, Morgane Griesbeck, Andrew Butler, Shiwei Zheng, Suzan Lazo, Laura Jardine, David Dixon, Emily Stephenson, Emil Nilsson, Ida Grundberg, David McDonald, Andrew Filby, Weibo Li, Philip L. De Jager, Orit Rozenblatt-Rosen, Andrew A. Lane, Muzlifah Haniffa, Aviv Regev, and Nir Hacohen. Single-cell RNAseq reveals new types of human blood dendritic cells, monocytes, and progenitors. Science (New York, N.Y.), 356(6335):eaah4573, April 2017. ISSN 1095-9203. doi: 10.1126/science.aah4573. [61] Binhuan Wang, Yilong Zhang, Will Wei Sun, and Yixin Fang. Sparse Convex Clustering. Journal of Computational and Graphical Statistics, 27(2):393–403, April 2018. ISSN 1061-8600. doi: 10.1080/10618600.2017.1377081. URL https://researchwith. njit.edu/en/publications/sparse-convex-clustering/. [62] Daniela M. Witten and Robert Tibshirani. A framework for feature selection in clustering. Journal of the American Statistical Association, 105(490):713–726, June 2010. ISSN 0162-1459. doi: 10.1198/jasa.2010.tm09415. URL https://www.ncbi.nlm.nih. gov/pmc/articles/PMC2930825/. [63] F. Alexander Wolf, Philipp Angerer, and Fabian J. Theis. SCANPY: large-scale single-cell gene expression data analysis. Genome Biology, 19(1):15, February 2018. ISSN 1474-760X. doi: 10.1186/s13059-017-1382-0. URL https://doi.org/10.1186/ s13059-017-1382-0. [64] Kai R. Wycik, Tiffany M. Tang, Tarek M. Zikry, and Genevera I. Allen. Cluster Analysis with Resampling for Validation and Exploration (CARVE), May 2026. URL http: //arxiv.org/abs/2606.00327. arXiv:2606.00327 [stat.ME]. [65] Eric P. Xing and Richard M. Karp. CLIFF: clustering of high-dimensional microarray data via iterative feature filtering using normalized cuts. Bioinformatics, 17(suppl_1): S306–S315, June 2001. ISSN 1367-4811, 1367-4803. doi: 10.1093/bioinformatics/ 17.suppl_1.S306. URL https://academic.oup.com/bioinformatics/article/17/ suppl_1/S306/262530. [66] Rui Xu and Donald C. Wunsch. Clustering Algorithms in Biomedical Research: A Review. IEEE Reviews in Biomedical Engineering, 3:120–154, 2010. ISSN 1937-3333, 1941-1189. doi: 10.1109/RBME.2010.2083647. URL http://ieeexplore.ieee.org/ document/5594620/. [67] Tianyi Yao and Genevera I. Allen. Feature Selection for Huge Data via Minipatch Learning. Technical Report arXiv:2010.08529, arXiv, February 2021. URL http:// arxiv.org/abs/2010.08529. arXiv:2010.08529 [cs, stat] type: article. 35

[68] Tianyi Yao, Daniel LeJeune, Hamid Javadi, Richard G. Baraniuk, and Genevera I. Allen. Minipatch Learning as Implicit Ridge-Like Regularization. In 2021 IEEE International Conference on Big Data and Smart Computing (BigComp), pages 65–68, January 2021. doi: 10.1109/BigComp51126.2021.00021. URL https://ieeexplore. ieee.org/abstract/document/9373110. ISSN: 2375-9356. [69] Bing Yin and Greg Hamerly. Hierarchical Stability-Based Model Selection for Clustering Algorithms. In 2009 International Conference on Machine Learning and Applications, pages 217–222, December 2009. doi: 10.1109/ICMLA.2009.64. URL https://ieeexplore.ieee.org/document/5381839. [70] Bin Yu and Karl Kumbier. Veridical data science. Proceedings of the National Academy of Sciences, 117(8):3920–3929, February 2020. doi: 10.1073/pnas.1901326117. URL https://www.pnas.org/doi/10.1073/pnas.1901326117. [71] Jesse M. Zhang, Govinda M. Kamath, and David N. Tse. Valid Post-clustering Differential Analysis for Single-Cell RNA-Seq. Cell Systems, 9(4):383–392.e6, October 2019. ISSN 2405-4712. doi: 10.1016/j.cels.2019.07.012. URL https://www.sciencedirect. com/science/article/pii/S2405471219302698. [72] Grace X. Y. Zheng, Jessica M. Terry, Phillip Belgrader, Paul Ryvkin, Zachary W. Bent, Ryan Wilson, Solongo B. Ziraldo, Tobias D. Wheeler, Geoff P. McDermott, Junjie Zhu, Mark T. Gregory, Joe Shuga, Luz Montesclaros, Jason G. Underwood, Donald A. Masquelier, Stefanie Y. Nishimura, Michael Schnall-Levin, Paul W. Wyatt, Christopher M. Hindson, Rajiv Bharadwaj, Alexander Wong, Kevin D. Ness, Lan W. Beppu, H. Joachim Deeg, Christopher McFarland, Keith R. Loeb, William J. Valente, Nolan G. Ericson, Emily A. Stevens, Jerald P. Radich, Tarjei S. Mikkelsen, Benjamin J. Hindson, and Jason H. Bielas. Massively parallel digital transcriptional profiling of single cells. Nature Communications, 8(1):14049, January 2017. ISSN 2041-1723. doi: 10.1038/ncomms14049. URL https://www.nature.com/articles/ncomms14049.

36