Cell neighborhood topology directs rare cell population identification

Domain: biology
Published: 2026-03-30
Editorial outcome: major_revision

Abstract
Advances in single-cell and spatial technologies have transformed the dissecting of cell composition and
tissue architecture in complex biological systems. However, identifying rare cells critical to disease
pathogenesis and biological processes remains difficult, as their low abundance often masks them among
dominant cell populations. We introduce RareQ, a fast and scalable framework for rare cell detection by
evaluating cliquishness of each cell’s k-nearest neighborhood using single-cell omics data. Extensive
benchmarking on diverse simulated and real datasets demonstrates that RareQ exceeds existing methods in
accuracy, sensitivity, and efficiency. RareQ also excels at identifying modality-specific and -shared rare cells.
Its versatility across various biological contexts enables the discovery of functionally distinct rare cells with
unique molecular signatures in physiological and pathological context. RareQ's application to spatial
transcriptomics data reveals anatomically distinct and clinically relevant rare cell populations. Together, RareQ
offers an efficient approach for rare cell discovery, enhancing insights into tissue organization and disease
mechanisms.

Introduction
Rapid development of single-cell and spatial omics technologies has revolutionized our ability to dissect
complex tissues at unprecedented resolution1-3. These technologies not only enable the identification of
diverse cell types and functional states, but also provide insights into spatial organization, cell–cell interactions,
and dynamic transitions within tissue contexts3. Despite substantial progress in detecting major cell
populations using clustering algorithms, accurately identifying and characterizing rare cell types remains a
significant challenge. Rare cell populations—often comprising less than 1% of all profiled cells4-7—frequently
play outsized roles in physiological and pathological processes8. For example, circulating tumor cells (CTCs)
and cancer stem cells (CSCs), though scarce, can initiate metastasis and sustain tumor growth9-11. In
neurological disorders, rare neural stem cells influence neurotoxic activation of reactive astrocytes12. Standard
clustering algorithms are biased toward main cell populations and often merge or discard small clusters13.
Accurate detection of rare cell types demands computational approaches that are robust to technical noise
and sensitive to weak signals of rare cells.

In response to these challenges, a number of methods have been proposed to improve rare cell
identification in single-cell RNA sequencing (scRNA-seq) data. Clustering-based approaches enhance
sensitivity to small populations. For instance, RaceID detects rare cell populations by identifying outliers with
significantly deviating gene expression profiles from cluster medoids14. GiniClust2 integrates Gini index and
Fano factor-based gene selection to identify both rare and common cell types15,16. CellSIUS adopts a two-
step framework involving initial clustering followed by subclustering using genes of bimodal expression within
pre-defined clusters17. Other strategies eschew clustering by treating rare cells as statistical outliers. FiRE
assigns each cell a “rareness” score based on its local density using a sketching-based hashing scheme to
efficiently identify sparsely populated regions18. Similarly, EDGE leverages ensemble learning and
dimensionality reduction to highlight consistent outliers across multiple low-dimensional embeddings19. As
datasets scale to millions of cells, computationally-efficient algorithms are introduced, such as GapClust, a
lightweight approach tailored for ultra-large scRNA-seq datasets, balancing speed and accuracy20. Meanwhile,
scCAD recursively decomposes large clusters based on internal heterogeneity, enabling the recovery of rare
subpopulations masked in initial partitions21. A critical limitation of current methods is their modality-specific
design—most are exclusively optimized for scRNA-seq and cannot generalize to other single-cell modalities
such as scATAC-seq or spatial transcriptomics. Methods like MarsGT aim to overcome this by integrating
RNA and ATAC modalities through a heterogeneous graph transformer22. While promising, their reliance on
paired multi-omics datasets limits applicability in practical settings where such data are unavailable. A
promising yet underexplored direction for rare cell discovery lies in leveraging the topological properties of
low-dimensional data embeddings.

Here, we introduce RareQ, a topology-aware framework that exploits a cell's local neighborhood
connectivity to sensitively and efficiently identify rare cell populations. We benchmark RareQ by comparing
with existing methods, across a wide range of simulated and real datasets, spanning diverse biological
contexts, technological platforms, and multi-modal data types. Its lightweight design makes it computationally
efficient and suitable for ultra-large datasets, where many existing methods fail due to memory or runtime
constraints. We also demonstrate that RareQ can integrate with multi-omics frameworks to detect rare cell
populations that might be missed when analyzing individual regulatory layers separately. Finally, RareQ
proves its value in five applications. With airway epithelial scRNA-seq data, it identifies rare cell populations
involved in cell cycle regulation, ciliogenesis, and response to interferon functions. In B lymphoma multiome
data, RareQ identifies rare myeloid cells such as cDC with enhanced PD-L1 expression. In Alzheimer’s
disease (AD) data, RareQ finds microglia and astrocyte populations resided in reactive states. In a spatial
transcriptomic data of mouse brain, RareQ uncovers small, anatomically significant subpopulations that were
overlooked by previous analyses and existing rare cell detection methods. In a human renal cell carcinoma
(RCC) spatial data, RareQ captured the tertiary lymphoid structures (TLS) and rare myeloid and proliferative
cell populations that were missed by existing algorithms. Taken together, these benchmarking results and
biological applications highlight the power and generalizability of RareQ in uncovering rare but biological
significant cell populations from increasingly large and complex single-cell and spatial datasets.

Results
Overflow of RareQ
The topological structures of single-cell data are typically represented as low-dimensional k-nearest neighbor
(kNN) graph which models cells as nodes and proximity relationships as edges. In this network, rare cells
tend to form small cliques: their neighbors are tightly interconnected, with strong internal links and few
connections to external cells, in contrast to abundant populations. To quantify this property, we defined a
metric termed as neighborhood connectivity (Q). Q measures the proportion of internal connections within a
cell’s neighborhood relative to the total number of connections involving that cell and its neighbors.
Conceptually, Q is similar to the clustering coefficient of “ small-world” network for measuring the cliquishness
of a friendship circle23, but here it is reformulated as a cell-specific topological measure that characterizes
each cell’s local neighborhood network. Accordingly, high-Q cells and their tightly knit neighborhoods can
serve as topological anchors around which homologous cells can be preferentially aggregated (Fig. 1a). This
metric can be extended to the cluster-level, denoted as Q , to evaluate cluster independence by calculating
the ratio of intra-cluster connections to the total number of connections.

As expected, simulated data analysis shows that rare cells tend to exhibit higher Q values compared to
cells from major populations, indicating that rare cells are more likely to form topologically distinct, internally
cohesive structures (Fig. 1b). At the population level, both rare and major clusters tend to have Q values
approaching 1, reflecting their strong internal connectivity. These observations suggest that the Q metric of
individual cells enables effective differentiation of rare but biologically relevant clusters, while its cluster-level
measurement guarantees the quality of identified clusters. Moreover, its ability to capture diverse topological
structures makes it particularly suitable for analyzing spatial transcriptomics data with complex spatial
architectures.

Based on the topological property of rare cell population, we developed a Q-based label propagation
strategy to infer both rare and major cell populations (Fig. 1c). The algorithm begins by assigning each cell to
its own initial cluster. During iterative propagation, each cell adopts the cluster label of the neighboring cell
with the highest Q value if its own Q is lower. In this way, high-Q rare cell cliques serve as stable topological
anchors that guide label diffusion outward across the kNN graph until convergence. To minimize
misclassification, we incorporate a kNN-voting mechanism that ensures cells connected to multiple clusters
are preferentially assigned to the one with greater connectivity. Clusters with average Q values exceeding a
predefined threshold are retained as high-confidence rare cell populations. For clusters with lower average Q
values, we apply a recursive merging strategy to identify major populations: if the merged cluster exhibits a
higher Q value than its constituent parts, the merge is accepted; otherwise, the clusters remain separated
(see Methods). We also performed ablation studies to isolate RareQ’s core components and evaluate their
individual contributions to rare-cell detection and overall clustering accuracy (Supplementary Note 1). The
results highlight Q-guided network propagation and post hoc merging as the central mechanisms underlying
RareQ, supporting its robust ability to identify rare cell populations across multiple data modalities and tissue
architectures.

Benchmarking RareQ in both simulated and real RNA-seq datasets
To assess the performance of RareQ, we compared it against seven existing rare cell detection methods (i.e.,
scCAD, CellSIUS, RaceID, FiRE, GiniClust2, EDGE, and GapClust) using three simulated experiments
encompassing 150 simulated scRNA-seq datasets. These datasets were generated by down-sampling
heterogeneous peripheral blood mononuclear cell (PBMC) data, each comprising 5,000 cells and representing
5 to 15 distinct cell types (Supplementary Table 1). The F score, which effectively captures the trade-off
between precision and recall, was used to evaluate the performance of different methods. Compared to other
methods, RareQ consistently surpassed them in detecting simulated rare cell populations (Fig. 2a). In Sim-
PBMC-1, RareQ outperformed the second-ranked method, scCAD, achieving a 3% improvement; In Sim-
PBMC-2, RareQ achieved a 7% improvement compared to the second-ranked method CellSIUS; In Sim-
PBMC-3, RareQ achieved 2%, 4%, 6% and 13% improvements compared to the second-best-performing
methods in identifying R1, R2, R4 and R5 rare cell types, respectively.

We further evaluated RareQ on the RNA modality of four human PBMC datasets (PBMC-bench-1, 2, 3
and 4) with known cell type labels (Fig. 2b and Supplementary Fig. 1). We defined rare cell clusters as those
comprising less than 1% of the total cell population. RareQ successfully identified 15 out of 17 rare clusters
across the four datasets, achieving the highest F scores among all methods (Fig. 2c). In contrast, the next-
ranked methods, scCAD and RaceID, found only 3 and 6 rare cell clusters, respectively. To confirm RareQ's
robust performance, we tested it on 20 additional independent scRNA-seq datasets from various tissues and
profiling platforms (Supplementary Table 2). These comparisons consistently showed that RareQ maintained
the highest F scores (median = 0.91), precision (median = 0.93), and recall (median = 0.91), significantly
outperforming other approaches (Fig. 2d and Supplementary Fig. 2). Notably, RareQ was able to identify most
rare cell types that other methods missed. RareQ identified 110 of 133 rare cell types, significantly
outperforming the second-ranked approach, scCAD, which found only 53 rare cell clusters (Fig. 2e). Across
various datasets, UMAP visualizations clearly demonstrated that rare cells identified by RareQ consistently
formed distinct clusters, separate from major cell populations (Supplementary Fig. 3-21). Furthermore, we
evaluated the global clustering performance of RareQ in identifying both major and rare cell populations using
the normalized mutual information (NMI) metric. RareQ consistently outperformed all competing methods
across simulated rare cell clusters. Notably, it achieved NMI improvements of approximately 17%, 7%, and
4% over the second-best method, CellSIUS, in three simulation experiments (Supplementary Fig. 22).
Additionally, we compared RareQ with competing methods across a wide range of parameter combinations,
and clustering resolutions. Across all settings, RareQ remained robust to parameter choices and achieved
superior performance relative to other methods (Supplementary Note 2-7).

To evaluate RareQ’s sensitivity and robustness to the number of differentially expressed (DE) genes, we
conducted experiments using an artificial scRNA-seq dataset and a Jurkat scRNA-seq dataset, each of which
comprised two cell types with the minor cell type representing approximately 1% of the total population. Genes
that were differentially up-regulated in the rare cell group served as DE genes (211 for the simulated dataset,
148 for the Jurkat dataset, see Methods). In each simulation, an equivalent number of non-DE genes was
replaced with randomly selected pre-identified DE genes. This process was repeated 10 times for each
number of DE genes, where the average F score was computed across iterations for comparison among
algorithms. It showed that RareQ precisely identified rare cell types, outperforming its competitors in both
datasets (Supplementary Fig. 23). While RaceID performed slightly better or comparably in the artificial
scRNA-seq dataset, its performance significantly declined in the Jurkat dataset. Furthermore, using the
aforementioned datasets, we confirmed that RareQ is highly stable, showing only minor performance
fluctuations near the boundaries of the tested parameter ranges (Supplementary Note 2, Supplementary Fig.
24-28). Batch-only data analysis demonstrated that applying batch correction is critical to prevent spurious
“rare populations” that are driven primarily by batch composition rather than biology (Supplementary Note 6).
Benchmarking on negative control data indicated that RareQ does not tend to invent rare clusters in relatively
homogeneous data when true rare populations are absent (Supplementary Note 7). Collectively, these
findings highlighted RareQ's robustness and its capacity for precise rare cell type identification, even when
differential expression signals among cell types were weak.

RareQ is fast and highly scalable
To evaluate the scalability of RareQ on large datasets, we applied it to a mouse brain scRNA-seq dataset
containing 1.47 million cells spanning 31 distinct cell types24. Among them, 15 cell types representing less
than 1% of the total population were defined as rare cell populations (Fig. 3a). Of the eight tested methods,
only RareQ, FiRE, and GapClust successfully processed the full dataset on a computer platform x86_64
architecture, 64 CPUs (Intel Xeon Gold 6226R, 2.90GHz ) and ~500 GB of memory; the others failed due to
memory error. RareQ successfully identified all rare cell clusters, achieving superior accuracy compared to
FiRE, whereas GapClust failed to detect these rare cells (Fig. 3b). Notably, RareQ exhibited the highest
computational efficiency, requiring the least memory (144 GB) and the shortest runtime (30 minutes),
outperforming FiRE (529 GB, 0.79 hours) and GapClust (206 GB, 1.6 hours) (Fig. 3c).

We further evaluated RareQ using a series of eight simulated scRNA-seq datasets ranging from 1,000 to
200,000 cells, each consisting of one rare cell cluster (10 cells) and one major population. RareQ processed
10,000 cells in just 27 seconds and scaled to 200,000 cells in 13 minutes, achieving performance comparable
to FiRE (13 seconds and 7 minutes) while significantly outperforming GapClust (180 seconds and 75 minutes)
(Fig. 3d, Supplementary Table 3). Other methods exhibited limited scalability: RaceID failed on datasets with
over 50,000 cells; GiniClust2 and EDGE consumed over 253 GB and 35 GB of memory in processing 100,000
cells, respectively, while both crashed on datasets exceeding 100,000 cells (Supplementary Table 4); When
processing 200,000 cells, scCAD and CellSIUS required over 10 hours and 32 hours, respectively. In contrast,
RareQ was the only method that successfully identified rare cells in datasets exceeding 100,000 cells.
Comprehensive comparisons across three evaluation metrics—F score, precision, and recall—further
confirmed RareQ’s high sensitivity and accuracy (Fig. 3e, Supplementary Fig. 29). Taken together, RareQ
exhibited exceptional speed and highest memory efficiency on large-scale datasets.

RareQ enables rare cell detection from multi-modal single-cell data
To evaluate the applicability of RareQ on multi-modal single-cell data, we applied RareQ to a CITE-seq dataset
of human bone marrow mononuclear cells, which includes gene expression profiles and measurements of 25
surface proteins across 30,672 cells spanning 27 cell types25. We aimed to evaluate RareQ's ability to detect
rare cell populations—which might be missed in single-modality analyses—by integrating it with the weighted
nearest neighbor (WNN) framework, a recently developed approach for combining multiple modalities
measured within the same cell25. WNN assigned cell-specific weights to each modality according to its relative
information content and generated a unified low-dimensional embedding that retains modality-specific signals
(Supplementary Fig. 30a). Notably, several T cell subtypes were indiscernible in the RNA modality but clearly
resolved in the protein modality. For example, CD8⁺ naive T cells, which showed high transcriptional similarity
to CD4⁺ T cells, were clearly separated in embedding space of the protein profiles.

Using curated cell type annotations as ground truth, we benchmarked RareQ across three settings: gene
expression alone, protein expression alone, and WNN-integrated embeddings. RareQ successfully identified
7 out of 9 rare cell types across all scenarios, however, some of them were modality-specific (Supplementary
Fig. 30b). Specifically, progenitor megakaryocytes (prog_MK) and hematopoietic stem cells (HSCs) were
detected in both gene expression data and WNN-integrated data, while CD56+ bright NK cells and regulatory
T cells (Tregs) were exclusively identified in the protein modality (Supplementary Fig. 30c). Compared to
existing algorithms, RareQ consistently achieved the highest F scores across all rare cell types, with the sole
exception of CD56+ bright NK cells, which only RaceID detected (Supplementary Fig. 30d). Collectively, these
results demonstrated that RareQ can be effectively applied to WNN-integrated multi-modal data to enhance
rare cell identification. However, some rare cell types were detectable only in individual modalities,
underscoring the importance of performing independent analyses prior to integration to maximize detection
sensitivity.

RareQ outperforms the multi-omics integration method MarsGT in rare cell inference
Recent efforts, such as MarsGT, have introduced deep learning frameworks designed to use both scRNA-seq
and scATAC-seq data to identify rare cell populations22. However, integrating data modalities that have
different dimensionalities and distributions can lead to alignment errors or overlooked modality-specific signals.
In contrast, RareQ flexibly supports both independent analysis of each modality and integrated analysis in
combination with the WNN framework. Firstly, we benchmarked RareQ against MarsGT using three simulated
experiments of PBMC datasets (Sim-PBMC 1–3), each comprising 50 paired scRNA-seq and scATAC-seq
profiles. With 150 simulation runs, RareQ's F scores were consistently and significantly higher than MarsGT's,
especially when utilizing WNN-integrated embeddings (Fig. 4a).

Next, we applied RareQ and MarsGT to the aforementioned four PBMC multi-omics datasets (PBMC-
bench 1-4) with paired RNA and ATAC data, each containing 2 to 6 rare cell types (<1% total population).
In PBMC-bench-1 and -2, RareQ accurately identified all rare cell types using both the RNA modality and the
WNN-integrated embedding, whereas MarsGT detected only one rare cluster in each dataset—MK/E prog in
PBMC-bench-1 and plasma cells in PBMC-bench-2 (Fig. 4b, Supplementary Fig. 31-32). In PBMC-bench-3
and -4, MarsGT failed to detect any rare cell types, while RareQ successfully identified one of the two rare
clusters in PBMC-bench-3 and all rare cell types in PBMC-bench-4 (Fig. 4b, Supplementary Fig. 33-34).
In addition to the annotated rare cell types, RareQ effectively discovered novel rare subpopulations that
possessed significant chromatin features distinct from canonical lineages. An example from PBMC-bench-4
showed RareQ identified two granulocyte/macrophage progenitor (G/M prog) clusters, two plasmacytoid
dendritic cell clusters (pDC), and two lymphoid progenitor (Lymph prog) cell clusters, none of which were
detected by MarsGT (Fig. 4c). Differential analysis of gene expression and chromatin accessibility confirmed
the distinctiveness of these subtypes, with particularly pronounced differences in their chromatin accessibility
profiles (Fig. 4d, Supplementary Fig. 35). Notably, in terms of computational efficiency, we observed that the
GPU-accelerated version of MarsGT accomplish analysis with more than 10 hours (10 epochs for training) on
the computer with 2 GPUs of NVIDIA RTX A6000 with 48 GB memory, while the CPU version required seven
days when using 128 CPUs (AMD EPYC 7543 32-Core Processor) to process the PBMC-bench-2 dataset
(~15,000 cells). In contrast, RareQ completed the same analysis in just one minute. RareQ demonstrated
consistently superior performance to MarsGT across ten independent paired multi-modal datasets, even when
MarsGT was run with its optimal hyperparameters and with WNN-derived embeddings (Fig. 4e,
Supplementary Fig. 36-46). This strong performance was also reflected in its robust global clustering
capabilities confirmed across a comprehensive collection of diverse datasets (Fig. 4f, Supplementary Fig. 47).
These results collectively highlight RareQ is capable of uncovering novel, biologically meaningful cell clusters.

RareQ enables comprehensive identification of rare airway epithelial cell types
While basal, secretory, and ciliated cells are the most common airway epithelial cells, rare epithelial subtypes
are increasingly recognized for their ability to sense environmental stimuli and initiate repair26,27. Their
dysregulation often contributes to disease development. To comprehensively identify rare cell populations
within the airway epithelium, we applied RareQ alongside competing methods to a scRNA-seq dataset from
mouse tracheal epithelium26 (Fig. 5a). Ten of RareQ's identified clusters, which had relatively small cell
populations, exhibited significantly higher Q values, allowing them to be clearly distinguished from more
abundant cell types (Fig. 5b). Among these, six of RareQ's identified clusters showed strong agreement with
the original study's annotations as indicated by high Jaccard similarity (Fig. 5c). Furthermore, the top five
upregulated genes in each RareQ-identified rare cluster matched the marker features of their corresponding
annotated cell types (Fig. 5d). Note that other competing methods failed to fully recover these rare populations.
For example, RaceID—the second-best-performing method—was unable to detect the Krt4/13⁺ cell
population (Fig. 5a).

In addition to the known rare cell clusters, RareQ also identified four novel rare clusters (C3, C4, C5, and
C10) that diverged from the original annotations. For example, cluster C3 showed decreased expression of
histone genes such as Hist1h2ap relative to canonical cycling basal cells (C6). Cluster C4 exhibited markedly
elevated expression of genes like Ccdc67 and Ccno compared to abundant cell cluster C11 (Fig. 5e). C4 also
strongly expressed genes related to chromosome condensation, spindle midzone assembly, and intraciliary
transport (Supplementary Fig. 48). Clusters C5 and C10, which resembled secretory clusters C12 and C13,
were marked by strong upregulation of interferon-stimulated genes such as Ifit3, Ifit3b, and Isg15, as well as
Foxj1, Lrrc23, and Tuba1a (Fig. 5e). Pathway-level analysis revealed that C5 upregulated genes associated
with interferon-α and interferon-β response pathways, while C10 showed higher enrichment for pathways
associated with intraciliary transport, axoneme assembly, and early endosome to Golgi transport
(Supplementary Fig. 48). In addition, we confirmed the presence of C3, C4, and C10 in specific regions of the
airway lumen across three independent human spatial transcriptomics datasets28,29, using the corresponding
RareQ-derived marker genes (Supplementary Fig. 49). Together, these results demonstrated that RareQ
enables comprehensive identification of rare airway epithelial cells and uncovers novel functional cell clusters.

RareQ identifies PD-L1 high dendritic cell subsets in B cell lymphoma
To evaluate the applicability of RareQ in cancer-related contexts, we analyzed 14,566 paired scRNA-seq and
scATAC-seq profiles obtained from intra-abdominal lymph node tumor tissue of a patient diagnosed with
diffuse small lymphocytic lymphoma. RareQ not only recapitulated the three rare populations (17, 18, and 16)
previously detected by MarsGT, i.e. B (13), CD4 T (11), and CD4 T (12) cells, but also discovered additional
rare populations, particularly within the macrophage lineage (Fig. 6a). These included clusters 5, 7, 9, 12, 15,
and 23, which displayed distinct molecular signatures (Fig. 6b).

Annotation with canonical marker genes classified clusters 9 and 23 as macrophages, clusters 5 and 7
as conventional dendritic cells (cDCs), and clusters 12 and 15 as monocytes (Fig. 6c). Clusters 5 and 7
exhibited pronounced antigen-presenting features, whereas clusters 9 and 23 were enriched for macrophage
activation signals (Fig. 6d). Importantly, cluster 5 demonstrated high expression of PD-L1 and activation of
the PD-L1/PD-1 checkpoint signaling pathway. Consistently, the transcription factors STAT1 and NFKB1—
key regulators of PD-L1 expression—showed markedly increased expression and chromatin accessibility at
their genomic loci (Fig. 6e-f and Supplementary Fig. 50). Given the critical role of PD-L1–high DCs in
dampening T cell activation and shaping responses to immune checkpoint blockade, precise identification of
this population may provide valuable insights for cancer immunotherapy. Additionally, RareQ identified rare
populations not detected by MarsGT, including fibroblasts (3), plasma cells (6), endothelial cells (8), and
plasmacytoid dendritic cells (pDCs, 11) (Supplementary Fig. 51).

RareQ uncovers rare microglial and astrocytic subsets in Alzheimer’s disease
To assess the utility of RareQ in neurodegenerative disease, we applied it to single-nucleus RNA-sequencing
(snRNA-seq) data from prefrontal cortex tissues of four healthy individuals and four patients with Alzheimer’s
disease30 (AD). RareQ resolved 25 distinct clusters, including five endothelial groups, seven astrocyte groups,
three microglial groups, four oligodendrocyte groups, and six oligodendrocyte precursor cell (OPC) groups
(Fig. 7a).

Within the microglial population, rare subsets 15 and 16 displayed distinct transcriptional programs (Fig.
7b). Cluster 15 was defined by elevated expression of inflammatory and phagocytosis-associated genes such
as IFI44L and NAMPT, whereas cluster 16 expressed neuronal connectivity–related genes including PCDH9
and IL1RAPL1. Abundance analysis revealed a strong enrichment of cluster 15 in AD patients (Fig. 7c,
Supplementary Table 6). This increase in C15 in AD was further confirmed by its signature gene expression
in an independent bulk RNA-seq cohort (Fig. 7d). Functionally, this subset showed enhanced activity in
phagocytosis, myeloid leukocyte activation, ion transport, and cytokine-mediated signaling pathways, all of
which are linked to the clearance of AD-associated pathological products (Fig. 7e). SCENIC31 analysis
revealed enhanced activity of transcription factors linked to pro-inflammatory regulation, including STAT1,
ELF1, ETV6, and FOXP1, as well as context-specific immune and phagocytic regulators such as RUNX2 (Fig.
7f).

Astrocyte heterogeneity was also evident across seven subpopulations (clusters 5, 7, 12, 18, 19, 20, and
24) (Fig. 7g). Clusters 5, 12, and 20 were preferentially enriched in AD patients, with cluster 5 consistently
detected across multiple individuals (Fig. 7h). All three expressed markers of reactive astrocytes; however,
cluster 12 exhibited a distinct molecular profile (Fig. 7i), suggesting alternative reactivation mechanisms under
pathological stress. ssGSEA analysis revealed that cluster 12 was significantly enriched for chaperone-
mediated autophagy and amyloid fibril formation pathways (Fig. 7j).

RareQ also identified rare endothelial, oligodendrocyte, and OPC subpopulations, each characterized by
distinct molecular signatures (Supplementary Fig. 52). In contrast, scCAD, CellSIUS, RaceID, GiniClust2,
FiRE, EDGE, and GapClust failed to detect microglial cluster 15 or astrocyte clusters 5 and 12 (Supplementary
Fig. 53).

RareQ identifies spatial architecture of rare cells from multi-platform spatial transcriptomics data
Advances in in situ hybridization (ISH) and in situ sequencing (ISS) technologies have enabled high-resolution
spatial transcriptomics through platforms such as Enhanced ELectric Fluorescence in situ Hybridization (EEL
FISH) and Xenium32,33. These platforms can map the spatial distribution of hundreds of RNA species at
subcellular resolution, allowing for the segmentation and identification of individual cells based on defined
gene panels. We evaluated RareQ’s capability to detect rare cell types and resolve fine-grained spatial
architectures using the spatial omics data. We first applied RareQ to an EEL FISH-derived spatial
transcriptomics dataset from a sagittal mouse brain section, comprising 127,591 cells and 440 genes
(Supplementary Fig. 54a). Notably, without relying on any external reference data, RareQ identified 13 rare
cell clusters (each comprising <1% of total cells) and revealed several small yet anatomically meaningful
populations corresponding to the rostral migratory stream (RMS), hippocampus, and the thin Purkinje layer of
the cerebellum32 (Supplementary Fig. 54b). The spatial expression of marker genes such as Sox11, Prox1,
and Gdf10 supported the anatomical identities of these cell clusters (Supplementary Fig. 54c). The side-by-
side comparison between RareQ-derived annotations and the original cluster labels was presented in
Supplementary Fig. 54d. By comparison, scCAD, CellSIUS, FiRE, and GapClust either did not detect rare cell
clusters or identified clusters that were diffusely distributed without correspondence to anatomical structures
(Supplementary Fig. 55).

RareQ also uncovered previously unreported rare-cell populations, exemplified by a choroid plexus
epithelial subgroup (C1) marked by high expression of Foxj1, Otx2, and Lratd2, and spatially embedded within
the major epithelial cluster (C11) (Supplementary Fig. 56a–c). The spatial positioning of C1, together with its
strong expression of cilium-associated genes such as Foxj1, suggests that this population corresponds to
ciliated epithelial cells of the choroid plexus. Because cilia play essential roles in cerebrospinal fluid flow and
maintenance of the ventricular microenvironment, identifying this rare ciliated-like epithelial population
provides a higher-resolution cellular framework for dissecting the structural and functional organization of the
ventricular system. We confirmed the presence of this cell type in an independent spatial dataset from early
mouse brain development34, which revealed a corresponding epithelial population expressing Foxj1 and Otx2
adjacent to the main Ttr-defined choroid plexus clusters (Supplementary Fig. 56d, e). Further validation using
a high-resolution choroid plexus scRNA-seq dataset demonstrated that epithelial populations could be
accurately mapped to their spatial counterparts in the EEL FISH data (Supplementary Fig. 56f, g). Among the
16 epithelial clusters identified by RareQ, several groups—including C3, C4, and C5—exhibited elevated
expression of Foxj1, Otx2, and Lratd2, mirroring the molecular signatures of the spatially defined C1
population (Supplementary Fig. 56h). Together, these findings strongly support that the rare populations
identified by RareQ within the spatial context represent authentic, specialized functional states rather than
technical artifacts.

We next applied RareQ to a second spatial transcriptomics dataset from a mouse brain coronal section,
generated using the Xenium platform. Compared to EEL FISH, Xenium provided lower inter-cell-type
contamination, with 162,033 cells profiled across 248 genes33. RareQ identified 199 cell clusters, of which
170 contained fewer than 1% of cells, reflecting its sensitivity to rare populations (Fig. 8a). We further
confirmed with CHOIR35-based testing that these clusters represent statistically distinct populations rather
than overclustering artifacts (Supplementary Fig. 57a). Furthermore, we used the silhouette coefficient—an
effective internal metric for assessing cluster count appropriateness and overall quality—to independently
evaluate RareQ’s clustering quality against Seurat (Supplementary Fig. 57b). It showed that RareQ's clusters
generally achieved a high-quality score than those generated by Seurat. The spatial organization of these
clusters revealed molecularly distinct anatomical regions, evidenced by differential gene expression when
compared to their most similar major clusters (Supplementary Fig. 58). In particular, RareQ delineated five
distinct cell clusters within the hippocampus, corresponding to CA1, CA2, CA3, and dentate gyrus (DG)
neurons. Notably, RareQ recovered the CA2 population, a region located between CA1 and CA3 that the
original data analysis had not uncovered. This cluster was characterized by the high expression of Stard536
and Plcxd3 (Fig. 8b). The CA1 and DG regions were confirmed by their specific markers, including Wfs1,
Fibcd1, Slc44a5, and Plekha2 (Fig. 8c). Interestingly, RareQ notably segmented the CA3 region into two
distinct subclusters (CA3_1 and CA3_2) that standard clustering algorithm failed to resolve (Supplementary
Fig. 59). Differential expression analysis revealed a gradient pattern of gene expression across the two
clusters. Specifically, Strip2 displayed a decreasing trend from CA3_1 to CA3_2, while Rspo2, Nwd2, and
Plcxd2 showed progressively increasing expression in the same direction (Fig. 8d). These findings highlighted
RareQ’s ability to uncover spatial heterogeneity that was not previously resolved.

To validate the existence of these two CA3 subpopulations, we further applied RareQ to an independent
hippocampal single-cell RNA-sequencing dataset comprising 12,418 cells37,38, which were partitioned into 76
clusters (Fig. 8e). The CA2 and CA3 populations were annotated using the established markers Cpne4,
Nectin3, and Plxcd3 (Fig. 8f). Within the CA3 subclusters, we identified distinct molecular signatures, including
Stard5, Strip2, Rspo2, and Nwd2, consistent with our spatial transcriptomic observations (Fig. 8g).
Comparative transcriptomic analysis revealed significant differences between CA2 and the three CA3
subclusters (Fig. 8h). ssGSEA-based pathway enrichment analysis further highlighted divergent activities in
amine binding, G-protein coupled purinergic receptor (GPCR) signaling, and phosphate ion binding across
these subpopulations (Fig. 8i). Moreover, previous evidence from RNA in situ hybridization and RNAscope
assays corroborated the molecular features of these two CA3 subtypes37,39-41. Besides, we validated our
findings using two independent mouse brain spatial transcriptomics datasets from 10x genomics datasets
(Supplementary Fig. 60).

Compared to RareQ, existing algorithms—including scCAD, CellSIUS, FiRE, and GapClust—exhibited
significant limitations in detecting rare spatial clusters. CellSIUS and GapClust failed to identify any rare spatial
clusters, while FiRE incorrectly detected rare populations with scattered, anatomically inconsistent
distributions (Supplementary Fig. 61). Although scCAD successfully identified the CA2 and CA3 rare cell types,
it missed the CA1 population and the subcompartments within CA3. Collectively, these results highlight RareQ
as a powerful and versatile tool for uncovering anatomical and molecular heterogeneity often overlooked by
existing methods.

RareQ identifies rare cells in spatial transcriptomics with paired protein data from renal cell carcinoma
To assess the efficacy of RareQ in analyzing multi-modal spatial data, we applied the algorithm to an in situ
gene and protein expression dataset derived from formalin-fixed paraffin-embedded (FFPE) human renal cell
carcinoma tissue available on the 10X Genomics website. This dataset encompassed 465,534 cells, profiled
for the expression of 405 genes and 27 proteins. Following integration using WNN analysis, RareQ delineated
59 distinct cell clusters, including 42 rare populations (Supplementary Fig. 62a). These clusters were distinctly
resolved in UMAP embeddings (Supplementary Fig. 62b). Annotation based on canonical markers identified
clusters corresponding to tumor cells, T cells, myeloid cells, B cells (including plasma cells), fibroblasts,
endothelial cells, and lymphatic endothelial cells (Supplementary Fig. 62c).

RareQ clustering revealed the presence of multiple tertiary lymphoid structures, comprising B cells, CD8⁺
T cells, CD4⁺ T cells, and endothelial cells, localized at the periphery of tumor regions (Supplementary Fig.
63a). Additionally, RareQ identified several rare myeloid subpopulations, including mast cells (clusters 2 and
4), conventional dendritic cells (cDCs; clusters 3 and 18), monocytes (cluster 6), and plasmacytoid dendritic
cells (pDCs; cluster 10) (Supplementary Fig. 63b). Notably, the cDC population in cluster 18 exhibited elevated
PD-L1 protein expression, consistent with our defined cDC subset in B-cell lymphoma. Furthermore, RareQ
detected several rare proliferative populations (Supplementary Fig. 63c), with marker-based annotation
identifying clusters 1 and 14 as proliferating tumor cells, cluster 5 as proliferating fibroblasts, cluster 7 as
proliferating CD4⁺ T cells, and cluster 12 as proliferating plasma cells (Supplementary Fig. 63d).

In the comparative analysis, scCAD and GapClust failed to identify these rare cell populations, while FiRE
detected only proliferative tumor cells. And CellSIUS identified certain cell populations that were poorly
resolved in UMAP embeddings (Supplementary Fig. 64). These findings highlight RareQ’s capability to identify
rare subpopulations in spatial transcriptomics with paired protein datasets.

Methods
Data preprocessing.
For scRNA-seq data, the raw count matrix was log2-normalized using Seurat R package42 (version 4.3.0).
Next, the top 2,000 highly variable features were selected and scaled for principal component analysis (PCA).
The k nearest neighbors (k.param, default 20) for each cell were determined using FindNeighbors function
(return.neighbor = TRUE) based on Euclidean distances computed in the embedded space of the top 50
principal components (PCs). Low-dimensional UMAP embeddings were also generated for visualization
based on the top 50 PCs.

For scATAC-seq data, the raw peak matrix was normalized using TF-IDF in Signac R package43 (version
1.9.0), and SVD was applied to derive a low-dimensional embedding. Like scRNA-seq data, the k.param
(default 20) nearest neighbors for each cell were determined based on Euclidean distances computed in the
SVD embedded space. Notably, the first SVD component was excluded from downstream analysis due to
high correlation with the number of fragments per cell43.

For spatial data from EEL FISH and Xenium platforms, Seurat R package was used for normalization,
scaling and PCA-based dimension reduction. Due to the low number of features (fewer than 500 genes) in
these data, we adjusted the scale factor to 80 and used all features as highly variable features in downstream
analysis. The k.param (default 20) nearest neighbors for each cell were determined based on Euclidean
distances in the low-dimensional PCA space.

For paired scRNA/ATAC-seq datasets, CITE-seq data, and spatial transcriptomics with paired protein
abundance data, standard analytical workflows were applied to each modality independently. WNN analysis
was then conducted to infer a WNN graph, denoting the most similar cells in the dataset based on a weighted
combination of protein and RNA similarities. Supervised principal component analysis (sPCA) was performed
based on the WNN graph, identifying the set of principal components that could transform the data in a single
modality to best capture the structure in a multimodal dataset. The k.param (default 20) nearest neighbors for
each cell were determined based on Euclidean distances in the low-dimensional sPCA space. The whole
process was accomplished using Seurat package.

For high-throughput single-cell data modalities beyond the workflow of Seurat, we prepared a helper
function to construct a pseudo-Seurat object based on the low-dimensional embeddings derived by
specialized tools. The k.param (default 20) nearest neighbors for each cell were determined based on
Euclidean distances in this low-dimensional embedding space.

RareQ then took the processed Seurat object above as input to infer rare cell types, and outputs the
predicted clusters. Thus, RareQ can be seamlessly incorporated into existing workflows by Seurat with a slight
parameter change, making it convenient for downstream analysis. Moreover, RareQ accommodates emerging
data modalities by constructing a pseudo-Seurat object from processed data.

RareQ algorithm
RareQ highlights rare cell signals via curated 𝑄 values, and further infers both major and rare cell populations
by propagating the labels of waypoint cells with high 𝑄 values to neighboring cells across the nearest
neighbor graph. Specifically, RareQ can be summarized into the following procedures after data
preprocessing.

1. Amplifying rare cell signals via 𝑄 metrics. Specifically, RareQ takes the processed Seurat object as input
and calculates the 𝑄 values for each cell with a neighborhood size 𝑘 of 6, where a small value
accommodates extremely rare cells. Let 𝐒={1,2,⋯,𝑘} be the set of 𝑘 cells including a cell and its 𝑘−1
nearest neighbors, its 𝑄 value can be calculated as:
∑ 𝐼
𝑄 = 𝑖,𝑗∈𝑺 𝑖,𝑗 (1)
𝑘×𝑘
1, An edge exists from 𝑖 to 𝑗
𝐼 ={ (2)
𝑖,𝑗 0, No edge exists from 𝑖 to 𝑗
𝐼 serves as an indicator, signifying the presence or absence of a connection from cell 𝑖 to 𝑗 in the kNN
𝑖,𝑗
graph. If 𝑖 =𝑗, 𝐼 =1. The denominator 𝑘×𝑘 measures the total number of directed edges from the 𝑘
𝑖,𝑗
cells to their 𝑘−1 nearest neighbors and themselves. 𝑄 ranges between 1⁄𝑘 and 1, a high value implies
dense internal connections within neighborhood.

2. Propagating cluster labels between neighboring cells by 𝑄 values. Initially, each cell is treated as its own
cluster. For each cell, it iteratively adopts the cluster label of its neighboring cell in 𝐒 with the largest 𝑄 value
if its own 𝑄 is lower. When a cell has multiple neighbors with the same highest Q value, RareQ breaks the
tie by selecting the neighbor with the smallest cell index. Once convergence is reached or the maximum
number of iterations is completed, the cluster labels of waypoint cells with local maxima of 𝑄 values are
propagated through the nearest neighbor graph, resulting in an initial set of clusters. Then the cluster label of
each cell is further refined by iteratively applying a majority-voted prediction among itself and its 𝑘−1
nearest neighbors.

3. Recursively merging cell clusters. Cell clusters with average 𝑄 values larger than 0.6 are retained as rare
cell populations, while the remaining clusters are subject to a merging procedure to infer main cell types.
Specifically, due to local variations of 𝑄 values, propagating labels of suboptimal waypoint cells might
introduce low-quality clusters, especially for main cell types. To merge these clusters, RareQ extends the 𝑄
metric to evaluate the internal connectivity within clusters, denoted as 𝑄 . For a given cluster C, its
𝑐
neighboring cluster C with the most connections from C is first identified. If the 𝑄 value of the combined
n 𝑐
cluster is greater than that of cluster 𝐶 and its neighboring cluster 𝐶 , and the ratio of the number of
n
connections from 𝐶 to 𝐶 to the total number of connections from cells in cluster 𝐶 exceeds 0.2, merging
n
is performed. Let 𝐂={1,2,⋯,𝑛} be a cluster with 𝑛 cells, its 𝑄 value can be calculated as:
𝑐
∑ 𝐼
𝑖,𝑗∈𝑪 𝑖,𝑗
𝑄 = (3)
𝑐 𝑛×𝑘.𝑝𝑎𝑟𝑎𝑚
1, An edge exists from 𝑖 to 𝑗
𝐼 ={ (4)
𝑖,𝑗 0, No edge exists from 𝑖 to 𝑗
𝑛, 𝑛 <20
𝑘.𝑝𝑎𝑟𝑎𝑚 ={ (5)
20, 𝑛 ≥20
Similar to 𝑄 metric, 𝐼 serves as an indicator, signifying the presence or absence of a connection from cell
𝑖,𝑗
𝑖 to 𝑗 in the kNN graph. If 𝑖 =𝑗 , 𝐼 =1. The denominator 𝑛×𝑘.𝑝𝑎𝑟𝑎𝑚 measures the total number of
𝑖,𝑗
directed edges from the 𝑛 cells to their 𝑘.𝑝𝑎𝑟𝑎𝑚−1 nearest neighbors and themselves. 𝑄 ranges
𝑐
between 1⁄𝑘.𝑝𝑎𝑟𝑎𝑚 and 1, a high value implies dense internal connections within a cluster.

Simulation of single-cell multi-omics data for benchmarking
Using paired human PBMC scRNA-seq and scATAC-seq data, we generated 150 simulated datasets following
a similar strategy adapted from Wang et al.22, with 50 replicates for each of the three simulation settings
(Supplementary Table 1).
Sim-PBMC 1 represented a baseline simulation setting. One cell type was randomly selected from the
two least abundant types in the original dataset to serve as the rare population, while the four most abundant
cell types were retained as major populations. Approximately 5,000 cells were sampled from these five cell
types; the rare population was fixed at 50 cells, and the remaining cells were allocated among the four major
populations according to their original proportions.
Sim-PBMC 2 increased the number of major populations to evaluate the method's performance in a more
complex background. Here, the nine most abundant cell types were retained as major populations, and one
was randomly selected from the two least abundant types as the rare population. As with Sim-PBMC 1, the
rare population was fixed at 50 cells, with the remaining cells sampled from the nine major populations based
on their original proportions, totaling approximately 5,000 cells.
Sim-PBMC 3 increased the number of both rare and major populations to assess the method's sensitivity
in detecting multiple rare cell types. In this setting, the ten most abundant cell types were retained as major
populations, and five additional cell types with frequencies below 1% in the original dataset were randomly
selected as rare populations. Approximately 5,000 cells were sampled from these 15 cell types according to
their original proportions.

In silico perturbation of scRNA-seq data for sensitivity evaluation
To analyze the sensitivity of RareQ to rare cell type identity, artificial scRNA-seq data was generated using
the splatter R package44. The following command was used to generate this data: splatSimulate(group.prob
= c(0.99, 0.01), method = ‘groups’, verbose = F, batchCells = 1500, de.prob = c(0.4, 0.4), out.prob = 0,
de.facLoc = 0.4, de.facScale = 0.8, nGenes = 5000, seed = 2024).
The dataset consists of 1,500 cells, each containing 5,000 genes. Of these cells, 1,487 represent the
major cell type, while 13 define the minor type. After data preprocessing, Wilcoxon’s rank sum test was used
to identify DE genes with an FDR cutoff of 0.05 and an inter-group absolute fold-change cutoff of 1.5. The
final 211 DE genes were preserved as a separate set. Meanwhile, another scRNA-seq dataset with the same
number of cells (1,500) and genes (5,000) was simulated with the above command but using de.prob = c(0,
0), which served as the non-differential gene set.
The subsampled Jurkat dataset consists of 1,556 cells, each containing 29,843 genes. Of these cells,
1,540 represent the major 293T cell type, while 16 define the minor Jurkat cell type. After data preprocessing,
Wilcoxon’s rank sum test was used to identify DE genes with an FDR cutoff of 0.05. To increase the number
of DE genes for analysis, the inter-group absolute fold-change cutoff was adjusted to 1. The final 148 DE
genes were preserved as a separate set. Additionally, 29,695 genes with a p-value exceeding 0.05 were
retained as a distinct set of non-differential genes. For ablation analysis, we also subsampled 293T cells from
major cell type and perturbed the expression of 10–100 genes by three-fold up regulation using the Jurkat
dataset.

Simulation of multi-scale scRNA-seq data for benchmarking
To analyze the scalability of RareQ to large data, a series of artificial scRNA-seq datasets were generated
using the splatter R package. Similar to the above command, the batchCells and group.prob parameters were
adjusted to generate datasets at different scales (1,000, 2,000, 5,000, 10,000, 20,000, 50,000, 100,000, and
200,000 cells) with only 10 rare cells in each dataset.

Statistics & reproducibility
In box plots, the horizontal line denotes the median value, the bottom and top of the box denote the lower and
upper quartiles that represent the 25th (Q1) and 75th percentiles (Q3) of values, respectively. The interquartile
range (IQR) is defined as the range between Q1 and Q3; whisker values are calculated as Q1 - 1.5 × IQR
and Q3 + 1.5 × IQR. Differentially expressed genes were identified using the FindMarkers function in the
Seurat R package (version 4.3.0), which employed a two-sided Wilcoxon rank sum test, with a false discovery
rate (FDR) cutoff of 0.05. Fold-change values were calculated based on the mean expression levels of each
gene between groups. P-values were adjusted using Bonferroni correction, accounting for the total number of
genes in the dataset. Differential abundance test between AD and healthy controls in microglia clusters were
applied using a negative binomial generalized linear model (NB-GLM) with quasi-likelihood (QL) estimation,
consistent with the approach implemented in the Milo45 algorithm.
No statistical method was used to pre-determine sample size. No data were excluded from the analysis,
all genes in datasets were used throughout all analyses. The selection of pre-identified differentially expressed
genes was randomized, all other experiments were not randomized. The investigators were blinded to
allocation during experiments and outcome assessment.

Comparison with other methods
Comparative analyses between RareQ and existing methods were conducted using default parameters. The
evaluated tools and their associated preprocessing pipelines are summarized as follows: (i) FiRE (v1.0.1,
GitHub), applied using default configurations. (ii) GapClust (v0.1.0, GitHub), data normalization was carried
out using the scran R package for most datasets. For the large-scale mouse brain dataset, normalization was
performed using Seurat to improve computational efficiency. (iii) CellSIUS (v1.0.0, GitHub), data were
normalized using scran, with an initial clustering step implemented via the quickCluster function. (iv) RaceID
(v0.3.5, CRAN), applied using default configurations. (v) GiniClust2 (v2.0, GitHub), applied with default
settings. (vi) EDGE (v1.0, GitHub), data preprocessing and clustering based on EDGE embeddings were
conducted using the Seurat package. (vii) scCAD (v1.0.0, GitHub), used with default settings as
recommended by the authors. (viii) MarsGT (v1.0.0, GitHub), for benchmarking on paired scRNA-seq and
scATAC-seq data from human PBMCs, the regulatory score matrix generated by MAESTRO46 was obtained
from Zenodo47.

Evaluation index
To evaluate the performance of existing rare cell identification tools, precision, recall, and F score metrics
were employed to quantify the accuracy in capturing each rare cell type. When evaluating the accuracy of
clustering-based algorithms like RaceID, GiniClust2 and CellSIUS that predict both major and rare cell types
simultaneously, we employed Normalized Mutual Information (NMI) to quantify the purity of the clustering
output. Moreover, Jaccard index was also applied to measure the similarity between annotated rare cell types
and predicted clusters. These metrics range from 0 to 1, where 1 provides the best performance. Before
evaluation, each predicted cluster was annotated by the majority type of its constituting cells. For each rare
cell type, the clusters with the same cell type annotation were considered as positive predictions.
Precision is defined as the proportion of correctly predicted rare cells relative to all rare cell predictions.
TP
Precision= (6)
TP+FP
Recall evaluates the proportion of correctly predicted rare cells from actual rare cells.
TP
Recall= (7)
TP+FN
F score is a weighted average of precision and recall:
Precision∗Recall
F score=2× (8)
Precision+Recall
TP (True Positive) denotes the number of cells that are correctly predicted as rare cells and are actually
rare cells. FP (False Positive) refers to the number of cells that are predicted as rare cells but are actually
common cells. FN (False Negative) represents the number of cells that are predicted as common cells but
are actually rare cells.
NMI measures the normalized dependency of the true labels on the predicted cluster, which was
calculated using the aricode R package (version 1.0.2).
Jaccard index quantifies the membership similarity between predicted rare cells and true rare cells for
each rare cell type. Let 𝐒={𝑠 ,𝑠 ,⋯,𝑠 } be the set of predicted clusters, and 𝐓={𝑡 ,𝑡 ,⋯,𝑡 } be the set
1 2 𝑆 1 2 𝑇
of true labels.
𝐒∩𝐓
Jaccard index= (9)
𝐒∪𝐓

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Although the authors present the RareQ method for detecting rare cell populations in single-cell and spatial data, this work still falls short of the conceptual, technical, and biological standards expected by Nature Communications. Crucially, the study did not include disease-related datasets (e.g., cancer or neurodegeneration), making it impossible to assess the generalizability of the method. In addition, the benchmark design is nearly identical to MarsGT, which is of questionable novelty and fairness. Finally, for publicly available spatial transcriptomics datasets, RareQ's annotations used for comparisons do not match those in the original study reports, suggesting the possibility of misannotation or overclustering. Taken together, these shortcomings prevent me from recommending the manuscript for publication.\n1. The manuscript initially specifies that rare cells are both functionally important and low in abundance, yet in the Results they are selected solely on abundance. This contradiction undermines the entire biological premise and renders subsequent claims unconvincing.\n2. All data presented in the manuscript are derived exclusively from healthy donor samples, with no performance evaluation on datasets related to cancer or neurodegenerative diseases. This absence of validation on disease-relevant samples significantly limits the generalizability and practical applicability of the proposed method.\n3. Ground-truth labels are taken at face value, and the study provides virtually no orthogonal evidence (e.g., immunostaining, lineage tracing, perturbation) to confirm the identity or significance of the purported rare populations. More regrettably, there is almost no in-depth functional investigation of the identified cell types, which substantially diminishes the scientific value and biological relevance of the study.\n4. The benchmark experiments presented in Figure 2 largely overlap with those used in the previous study (MarsGT). However, MarsGT is included only as a baseline for the multi-omics integration task and is notably absent from the main benchmark and other experimental comparisons. Given that the benchmark settings in this manuscript closely align with those in MarsGT, excluding MarsGT from broader evaluations is unjustified and undermines the credibility and completeness of the performance assessment.\n5. The separation of the two pDC subpopulations in Figure 4 relies on only subtle epigenomic differences, which are insufficient to justify defining them as distinct subtypes. Without clear evidence of functional divergence, the claim that these represent true subpopulations is unsubstantiated and premature.\n6. The biological interpretation presented is inappropriate. In Figure 6, cluster C10 is enriched for the pathway “protein localization to the presynaptic region.” Discussing the “presynaptic region” in the context of airway epithelial cells is highly unusual and seemingly irrelevant, which significantly undermines the reliability of the conclusions.\n7. The analysis of the mouse brain (EEL FISH) dataset in Figure 7a–c presents several serious issues:\n(1) RareQ’s clustering results merely reproduce the specific structures highlighted in the original study, rather than uncovering novel biological insights. Furthermore, the annotation is less detailed—for instance, hippocampal subregions such as CA1/2 and CA3 are not distinguished, as they are in the original reference.\n(2) Additionally, no baseline methods are included for comparison on this dataset, making it difficult to evaluate the distinct advantages or effectiveness of RareQ in this context.\n8. The discovery of rare cell types in Figure 7d appears to be driven largely by over-clustering. In the Xenium dataset, 162,000 cells were divided into 199 clusters—most containing fewer than 1% of cells. This suggests that the identification of rare populations stems more from arbitrary clustering granularity than from any unique strength of the method. Standard clustering tools could yield similar results, and the approach is highly sensitive to minor cellular heterogeneity, casting serious doubt on the biological validity of these rare cell types.\n9. The description of the methodology on page 13, lines 468–477, is notably unclear and problematic.\n(1)The parameter k mentioned here should be clearly distinguished from the k value defined in the Data Preprocessing section to avoid confusion.\n(2)According to the manuscript’s logic, since all cells in the set S are connected as nearest neighbors, this would imply that every cell is connected by edges, resulting in all calculated Q values being equal to 1.\n10. Accurate identification of rare cells should be based on a correct overall clustering performance. However, the manuscript lacks evaluation metrics assessing the quality of the global clustering results.\n11. Given the discussion on memory usage and runtime, it is necessary for the authors to specify the exact GPU and CPU models used in their experiments.\n12. Since the study only uses spatial transcriptomics data for testing, the umbrella term “spatial omics” is misleading and should not be used. The authors should add at least one additional spatial modality, for example, spatial ATAC-seq, to substantiate multi-omic applicability.\n13. In line 468, the parameter k is set directly to 6. Since k is a critical parameter, its impact on the experimental results should be explicitly discussed, with a dedicated analysis of how varying k affects performance.\n14. The absence of example code in the provided code repository hinders reproducibility and makes it difficult for users to effectively apply the method.", "reviewer_2": "RareQ is based on the idea of measuring fractional local connectivity (the so-called Q values in a k nearest neighbour sense) and associating rare groups to isolated clusters with high proportion of internal connections. This is a simple concept but practically sensible.\nClustering proceeds first by iterative assimilation through cluster label propagation and majority voting. This is then followed by a merging process.\nThe authors should comment on the sensitivity to key choices of parameters and how the values used in the analyses were chosen, specifically:\n1. The average Q value threshold of 0.6 used to retain rare cell populations.\n2. The ratio of connections from C to Cn to the total number of connections from cells in cluster C exceeding 0.2.\n3. The k = 20 threshold.\nThe authors provide source code on Github and analysis data and scripts in Zenodo. However, I could not identify the scripts for the simulation analysis. It would be beneficial if the full simulation analysis was made available to give complete confidence in the results presented.\nIt would be beneficial for the authors to provide tables precisely describing the number of cell populations and numbers in their simulation scenarios. This is currently a little unclear from the descriptions in the manuscript.\nOverall, the results appear promising but the reproducibility aspect of the work needs to be significantly stronger.", "reviewer_3": "Identifying rare cell types from single-cell omics data has been a challenging task. This paper introduces a simplisitic approach by propagating cluster labels on a knn graph according to the clustering coefficient, which seems to achieve significantly improved performance on the presented results. I find the description some of the presented evaluation metrics unclear, and it is hard to tell the general applicability and performance of the method due to a lack of 1) evaluations considering all cell types, and 2) a benchmark with single-cell data from diverse systems. Please see my detailed comments below.\nSpecific points:\n1.The “Q value” is the same as the very widely used “clustering coefficient”. Proper citation of the original paper (Watts, Duncan J., and Steven H. Strogatz. \"Collective dynamics of ‘small-world’networks.\" nature 393.6684 (1998): 440-442.) should be added. Throughout the article, it should not be described as a novel metric. Rather, it is more appropriate to put this as an application of a generally applicable classical index to the problem of single-cell data analysis.\n2.It is described in Methods-RareQ algorithm-section 2, that the label propagation is stopped until convergence or a maximum number of iterations is reached. Since there are thousands of nodes in typical single-cell/spatial data but only 6 possible distinct Q values, I guess it is quite common for a node to have several highest Q value neighbors. I am therefore wondering that in this case, how do you break the tie? If a non-determinstic criteria is used, so that the convergence is then not guaranteed, is the result sensitive to the number of maximum iterations? If a deterministic criteria is used (e.g., always choose the cell with smaller index), is the result sensitve to criteria (in the example, will the result change much with perturbed indices)?\n3.In the description of Evaluation index in the Methods section, the F1 score evaluates a binary classification task of whether a cell is a rare cell type or common cell type, according to lines 615 to 617. But the results (e.g., Fig.2) report cell type specific F1 scores. How are these computed exactly?\n4.On line 605, it is described that “the clusters with the same cell type annotation were considered as positive predictions”. In the evaluation, do you consider the situation where cells from the same rare cell type fall in different clusters? When there are multiple clusters whose major cell type are the same, are they all considered positive predictions? Since each cluster is annotated by the known cell type of majority of its consituting cells, in the extreme case, a clustering where each cell itself forms a cluster will all be considered positive predictions. And in a less extreme case, a high performance can be achieve under this metric by producing many small clusters.\n5.Related to the previous point, it is important to report overall clustering performance as well, such as NMI over all cell types. In the current manuscript, NMI is reported on rare cell types only (Supp. Fig. 1).\n6.Can the authors investigate more into whether the improved performance comes from the RareQ algorithm or due to better preprocessing of the data which led to higher quality distance metric among the cells? An added comparison to other methods using the the same preprocessing pipeline that RareQ used would be helpful to clarify this point.\n7.For the spatial application, does the knn graph only rely on gene expression similarity? It would be interesting to explore the usage of both gene expression similarity and spatial proximity to reveal “spatially rare” cell types. That is, a spatially distinct region of a small number of cells from a type that is considered common from the single-cell perspective.\n8.I appreciate the in-depth studies presented in the manuscript each focusing on one dataset (a total of three single-cell systems), for example, the hundreds of single-cell data simulated from the PBMC dataset. However, broad benchmark is also important to demonstrate the general applicability of RareQ. I suggest that the authors add a more comprehensive benchmark that includes real scRNA-seq datasets from diverse systems, using the numerous expert annotated single-cell data. Alternatively, at least including the datasets from the compared methods can be used to enrich the evaluation."}, "author_response": "We thank the reviewer for these critical comments. RareQ's conceptual basis is the Q value, a metric specifically designed to quantify the cliquishness (or interconnectedness) within a cell's k-Nearest Neighbors (kNN) network. Reviewer #3 highlights that this underlying concept is similar to the clustering coefficient used in a small-world network model (Watts & Strogatz, Nature, 1998; Kleinberg, Nature, 2000). Crucially, the novelty of our work is the adaptation of this classical network principle into a quantitative framework specifically for rare cell detection, not the benchmark design itself—since all methods are evaluated using the same benchmarking strategies and metrics for fairness.\nIn this revision, we have substantially expanded our evaluations to include a wide range of datasets, particularly disease-related datasets—including paired single-cell RNA and ATAC sequencing dataset from B cell lymphoma and single-nucleus RNA-seq dataset from Alzheimer’s disease, as well as spatial omics datasets from renal cell carcinoma—to rigorously assess the generalizability of RareQ. We have carefully addressed all concerns, including the fairness of comparisons with MarsGT, annotation consistency, and potential over-clustering in public spatial datasets, as detailed in the point-by-point responses below.\nWe believe these comprehensive revisions address the reviewer’s concerns and strengthen both the conceptual and biological contributions of this study. We respectfully ask the reviewer to reconsider and re-evaluate our revised manuscript.\n1. The manuscript initially specifies that rare cells are both functionally important and low in abundance, yet in the Results they are selected solely on abundance. This contradiction undermines the entire biological premise and renders subsequent claims unconvincing.\nResponse: We appreciate this insightful comment, which helps us clarify the biological rationale and scope of our study on rare cells. We recognize that our initial use of 'functionally important' may have created the impression that we experimentally validated the specific role of rare cells, which is not the primary focus of this computational study. To address this, we have refined our terminology: rare cells are now defined as subpopulations characterized by both low abundance and biological significance. The biological significance of these cells is evidenced by their distinct molecular profiles (e.g., unique transcriptomic, proteomic, or epigenomic signatures) that robustly differentiate them from dominant cell types. These unique profiles often suggest specialized roles, such as developmental intermediates, disease-driving subclones, or niche-specific populations. It is important to note that our method, RareQ, does not select cells based on a simple abundance threshold. Instead, it identifies cells based on their neighborhood relationships in low-dimensional space. This computational strategy ensures that the detected cells possess molecularly distinct gene expression profiles or occupy unique spatial niches, irrespective of their quantity.\nTo substantiate this biological relevance, we present extensive molecular evidence:\n Molecular Signatures: We conducted differential gene expression analyses and performed Gene Set Enrichment Analysis (GSEA) to link the rare cell signatures to biologically relevant processes and pathways.\n Disease Relevance: We included analyses of disease-related datasets (e.g., tumor and neurodegenerative contexts) that demonstrate the broad significance of the identified rare cells across diverse biological systems.\n Independent Benchmarking: We employed independent benchmarking methods, such as CHOIR, to verify that the detected rare cells are not statistical artifacts but exhibit molecular hallmarks consistent with specialized roles.\nTaken together, these analyses demonstrate that RareQ detects rare cells with molecular features indicative of biological significance, thereby providing a crucial foundation for subsequent functional and mechanistic studies.\n2. All data presented in the manuscript are derived exclusively from healthy donor samples, with no performance evaluation on datasets related to cancer or neurodegenerative diseases. This absence of validation on disease-relevant samples significantly limits the generalizability and practical applicability of the proposed method.\nResponse: Thank you for highlighting the need to evaluate RareQ on disease-relevant datasets to enhance its generalizability and practical applicability. In response, we have extended our analysis to include three diverse pathological datasets, demonstrating RareQ’s robust performance in detecting biologically significant rare cell populations in cancer and neurodegenerative contexts. These results are now incorporated into the revised manuscript.\n(1) Lymphoma Dataset: We applied RareQ to a matched single-cell RNA and ATAC sequencing dataset from a flash-frozen intra-abdominal lymph node tumor of a patient with diffuse small lymphocytic lymphoma (also used in MarsGT’s paper). RareQ identified rare cell populations with distinct marker genes—including plasmacytoid dendritic cells (pDCs), fibroblasts, endothelial cells, and plasma cells (Fig. R1), substantially exceeding the performance of MarsGT (Fig. R2a). We focused on rare populations within the myeloid lineage: clusters 5 and 7 corresponded to cDC subtypes, clusters 12 and 15 to monocytes, and clusters 9 and 23 to macrophages (Fig. R2b, c). Notably, cluster 5 displayed elevated PD-L1 expression and PD-1 checkpoint signaling activity, a feature previously linked to tumor immune suppression and checkpoint therapy response (Fig. R2d). Furthermore, both PD-L1 and its upstream transcriptional regulators STAT1 and NFKB1 exhibited high expression and chromatin accessibility (Fig. R2e, f and Fig. R3), supporting their key regulatory role driving the transition of cDCs toward an immunotolerant phenotype.\nFig. R1. Dot plot illustrates the expression values and proportions of marker genes in clusters predicted by RareQ in B lymphoma dataset.\nFig. R2 RareQ identifies rare myeloid cell types in B lymphoma. a, UMAPs visualizing cell clusters predicted by RareQ and MarsGT, with the common rare cell types identified by both highlighted. b, Dot plot demonstrates the expression value and proportion of the top differentially expressed genes between the myeloid cell clusters predicted by RareQ. c, Violin plot shows the expression of marker genes within the myeloid clusters for annotation. d, Gene signature scores of antigen presentation, macrophage activation, and PD-L1 pathways. P-values were calculated using ANOVA. Coverage plots for gene (e) STAT1 and (f) NFKB1.\nFig. R3 Coverage plots for gene PD-L1 on B lymphoma data.\n(2) Alzheimer’s Disease Dataset: Using single-nucleus RNA-sequencing data from cortical samples of non-cognitive impairment (NCI) controls and Alzheimer’s disease (AD) patients, RareQ identified multiple rare microglial and astrocytic subpopulations (Fig. R4a). Among microglia (marker genes CSF1R, SPP1), cluster 15 showed high expression of IFI44L and NAMPT, and was predominantly enriched in AD patients (Fig. R4b, c). Functional enrichment revealed strong activity of phagocytosis, myeloid leukocyte activation, and cytokine-mediated signaling pathways (Fig. R4d), suggesting an activated state with enhanced phagocytic capacity relevant to clearance of AD pathological products. Transcription factor analysis further showed enrichment of STAT1, RUNX2, and ELF1 (Fig. R4e), consistent with an activated microglial phenotype. In addition, RareQ detected several rare astrocyte subpopulations (clusters 5, 12, 20) that were disproportionately enriched in AD and displayed highly reactive molecular features (Fig. R4f-i), indicating pathological activation. In contrast, alternative algorithms (scCAD, CellSIUS, RaceID, GiniClust2, FiRE, EDGE, and GapClust) failed to accurately identify these biologically significant populations (Fig. R5), underscoring RareQ’s superior sensitivity.\nFig. R4 RareQ identifies rare activated microglia and astrocyte populations in Alzheimer’s Disease (AD). a, UMAP projections visualizing the major cell types and cell clusters predicted by RareQ. b, Dot plot illustrating the expression levels and proportions of top differentially expressed genes (DEGs) across microglial clusters predicted by RareQ. c, Bar plot showing the proportions of cells derived from AD patients and healthy control (CTRL) individuals within each microglial cluster. d, Bar plot displaying the enriched pathways of microglial clusters 15 and 16. e, Heatmap representing the inferred regulon activities across microglial clusters. f, Same as (b), but for astrocyte clusters. g, Same as (c), but for astrocyte clusters. h, Violin plot showing the expression of marker genes (reactive and homeostatic) within astrocyte clusters. i, Gene signature scores for the chaperone-mediated autophagy and regulation of amyloid fibril formation pathways. P-values were calculated using ANOVA.\nFig. R5 UMAP projections of the AD dataset colored by rare cell clusters predicted by the specified methods. Rare cell clusters (< 1% population) predicted by RareQ are highlighted for comparison.\n(3) Renal Cell Carcinoma Dataset: In FFPE human renal cell carcinoma samples with in situ gene and protein (Antibody Derived Tag, ADT) expression data, RareQ discovered several spatially distinct rare cell types, including immune-cell aggregates forming tertiary lymphoid structures, myeloid subpopulations, and diverse proliferative subtypes (Fig. R6, 7). Comparative analyses again confirmed RareQ’s advantages in detecting rare, spatially organized cell types within disease-associated multi-omics data (Fig. R8).\nThese findings demonstrate RareQ’s ability to robustly identify rare cell populations with biological and clinical relevance across diverse pathological contexts. We have added these results to the Results and Discussion sections, strengthening the manuscript’s claims regarding RareQ’s applicability and impact.\nFig. R6 Clusters predicted by RareQ in a multiome (RNA and ADT) spatial data from human renal cell carcinoma (RCC) by the Xenium platform. a, Spatial maps of RCC data, where every dot is a single cell colored by clusters predicted by RareQ, with selected rare cell clusters highlighted. b, WNN integration-based UMAP plot of RCC data colored by RareQ clusters. c, UMAP projections illustrating the expression levels of cell markers for annotation.\nFig. R7 RareQ identifies rare myeloid and proliferative cell types in multiome spatial RCC data. a, Zoomed view of spatial maps of RCC data showing tertiary lymphoid structures (TLS) colored by RareQ clusters and surface protein markers. b, Dot plot illustrating the expression values and proportions of marker genes and proteins in rare clusters in the myeloid compartment. c, UMAP plot of rare cell clusters with high expression of Ki-67 protein highlighted. d, Dot plot illustrating the expression values and proportions of marker genes and proteins in highlighted rare clusters in (c).\nFig. R8. UMAP projects of multiome spatial RCC data colored by rare cell types predicted by CellSIUS, FiRE and GapClust\n3. Ground-truth labels are taken at face value, and the study provides virtually no orthogonal evidence (e.g., immunostaining, lineage tracing, perturbation) to confirm the identity or significance of the purported rare populations. More regrettably, there is almost no in-depth functional investigation of the identified cell types, which substantially diminishes the scientific value and biological relevance of the study.\nResponse: We thank the reviewer for emphasizing the importance of validating the identified rare cell populations. We agree that orthogonal experimental approaches—such as immunostaining, lineage tracing, or perturbation studies—would provide definitive confirmation of their identities and biological significance. However, as such experiments require specialized models, resources, and conditions beyond the scope of this computational study. We now provide a more comprehensive analysis of molecular signatures, enriched biological pathways, and corroborating evidence from published experimental studies to support the identities and relevance of the rare populations uniquely detected by RareQ. For example, for the two rare CA3 subpopulations identified in Xenium mouse brain data, we reference independent studies employing in situ hybridization, RNAscope, and riboprobe-based methods, all of which confirm the specificity of the marker genes identified by RareQ (Fig. R9a-f). Moreover, we applied RareQ to analyze an independent scRNA-seq data (GSE261148) from mouse hippocampus, and successfully reidentified the CA3 subtypes with distinct expression profiles (Fig. R10a-e). Notably, markers (Stard5, Strip2, Rspo2, and Nwd2) discovered in the Xenium mouse brain data showed similar expression patterns within the CA3 subtypes in this scRNA-seq dataset (Fig. R10c).\nThese additional analyses strengthen the scientific foundation of our study and underscore RareQ's utility in uncovering rare populations that warrant future experimental investigation.\nFig. R9 Representative image of maker genes Rspo2, Nwd2 and Plcxd2 in published studies. a-d, Representative images of Rspo2 in situ hybridization in mouse hippocampus. e, Mouse hippocampus stained with a Nwd2 anti-sense riboprobe. f, Mouse coronal brain sections from the hippocampus probed for Plcxd2 (magenta) using RNAScope. The original publications were listed on top of plots.\nFig. R10 a, UMAP projects of cells from the mouse hippocampus colored by RareQ clusters, with CA2 and CA3 subtypes highlighted. b, Expression level of marker genes specific to CA2 and CA3 compartments. c, Violin plots showing the expression level of marker genes identified in Xenium mouse brain data across the CA2 and CA3 subtypes. d, Heatmap displaying the expression level of the top 5 DEGs between CA2 and CA3 subtypes. e, Violin plots visualizing the pathway enrichment score of selected pathways across the CA2 and CA3 subtypes. P-values were calculated using Wilcoxon test.\n4. The benchmark experiments presented in Figure 2 largely overlap with those used in the previous study (MarsGT). However, MarsGT is included only as a baseline for the multi-omics integration task and is notably absent from the main benchmark and other experimental comparisons. Given that the benchmark settings in this manuscript closely align with those in MarsGT, excluding MarsGT from broader evaluations is unjustified and undermines the credibility and completeness of the performance assessment.\nResponse: MarsGT was specifically designed for multi-omics datasets and, as noted in its original publication, cannot be applied to single-omics data. This inherent limitation precluded its inclusion in our initial single-omics benchmarks. We have previously benchmarked RareQ against MarsGT in paired PBMC simulated and real datasets (Fig. R11a, b). To provide a more comprehensive evaluation, and ensure a fair and rigorous comparison, we have now incorporated 10 additional multi-omics datasets and performed benchmarking using multiple metrics. The results show that RareQ consistently outperforms MarsGT across various metrics, including F score, precision, recall, and Normalized Mutual Information (NMI) (Fig. R11c, d). These expanded analyses underscore RareQ’s superior robustness and effectiveness in multi-omics rare cell detection. We have presented these complete benchmarking results in the revised manuscript (Fig. 4).\nFig. R11 Benchmarking of RareQ on multi-omics single-cell data. a, Benchmarking RareQ against MarsGT in the paired Sim-PBMC 1, 2, and 3 datasets via F score. RareQ was applied to each modality independently (RareQ_RNA for RNA data, RareQ_ATAC for ATAC data, RareQ_WNN for WNN-integrated data), while MarsGT integrates both modalities. b, Comparative results of RareQ against MarsGT on the four human PBMC paired scRNA-seq datasets (PBMC-bench-1, 2, 3 and 4) with scATAC-seq data via F scores. Benchmarking RareQ against MarsGT using 10 paired scRNA-seq and scATAC-seq datasets on (c) rare cell detection via F score, precision and recall, and (d) global clustering via NMI.\n5. The separation of the two pDC subpopulations in Figure 4 relies on only subtle epigenomic differences, which are insufficient to justify defining them as distinct subtypes. Without clear evidence of functional divergence, the claim that these represent true subpopulations is unsubstantiated and premature.\nResponse: We would like to clarify that the original Figure 4 only showed the top five loci with epigenomic differences. In the revised manuscript, we have now integrated transcriptomic, epigenomic, and regulatory analysis data to provide robust evidence for their distinct identities. Beyond the initially highlighted loci, our analysis identified 3618 differentially accessible regions (p < 0.01, absolute log fold change (|log FC|) > 1) (Fig. R12a, b). The scRNA-seq analysis further revealed 48 differentially expressed genes (p < 0.01, |log FC| > 1), including IL3RA, KLF4, KLF6, BACH2, FYN, and NRP1, that closely associated with pDC differentiation, function, or immune regulation (Fig. R12a).\nIn addition, SCENIC analysis identified distinct regulons, such as IKZF1, RFX5, REL, and FOXO3, which have clear associations with pDCs differentiation and function (Fig. R12c). Together, these detailed transcriptomic, epigenomic, and regulatory findings—now presented in Figure 4 and Supplementary Figure 35—provide strong evidence supporting the classification of these as genuinely distinct pDC subtypes.\nFig. R12 Molecular comparison of progenitor cell types in the PBMC-bench-4 dataset. a, Volcano plots of differentially expressed genes (RNA) and peaks (ATAC) between pDC_2 and pDC_1. b, Heatmap showing the top 100 differential peaks between pDC_2 and pDC_1. c, Heatmap visualizing the regulon activities of progenitor cell types.\n6. The biological interpretation presented is inappropriate. In Figure 6, cluster C10 is enriched for the pathway “protein localization to the presynaptic region.” Discussing the “presynaptic region” in the context of airway epithelial cells is highly unusual and seemingly irrelevant, which significantly undermines the reliability of the conclusions.\nResponse: We thank the reviewer for highlighting this important point. We agree that this pathway is biologically irrelevant to airway epithelial cells and may have caused confusion. This artifact likely arose because some genes participate in multiple pathways, which can lead to the statistical enrichment of processes not directly related to the system under study.\nTo improve the biological interpretability of our results, we manually filtered the initial enrichment outputs to exclude pathways that, although statistically significant, are biologically implausible in the context of airway epithelial cells. This included removal of the “protein localization to the presynaptic region” pathway (Fig. R13).\nFig. R13. Heatmaps displaying the top-enriched pathways of clusters within ciliated and secretory compartments, respectively\n7. The analysis of the mouse brain (EEL FISH) dataset in Figure 7a–c presents several serious issues:\n(1) RareQ’s clustering results merely reproduce the specific structures highlighted in the original study, rather than uncovering novel biological insights. Furthermore, the annotation is less detailed—for instance, hippocampal subregions such as CA1/2 and CA3 are not distinguished, as they are in the original reference.\nResponse: We acknowledge that in this dataset RareQ did not generate new biological discoveries. The original study reported 188 highly fine-grained subpopulations, distinguishing CA1/2 from CA3 within the hippocampus, which RareQ did not separate. We speculate that the authors may have used additional information to aid in their classification, as we were unable to reproduce these distinctions even when using more sensitive parameters to re-cluster CA subtypes. We've repeatedly tried to reach the authors for clarification, but haven't received a response.\nUpon closer examination, we found that the reported differences between the original CA3 and CA1/2 populations were minimal, with Cpne4 being the only gene showing a statistically significant difference (Fig. R14a, b). Moreover, benchmarking the original clustering with CHOIR (Cluster Hierarchy Optimization by Iterative Random forests) (https://doi.org/10.1038/s41588-025-02148-8), a recently published statistical framework specifically designed to guard against both over- and under-clustering in single-cell data, indicated that many of these fine-grained clusters likely reflected over-clustering rather than biologically distinct populations (Fig. R14c).\nDespite this, RareQ still outperformed existing rare cell detection methods on this dataset, demonstrating clear technical advantages (Fig. R15). Since no novel biological insights were derived here, we have moved these results to the supplementary file (now Supplementary Figure 52).\nFig. R14. Molecular comparison of CA1 and CA3 clusters in EEL FISH mouse brain spatial data. a. Volcano plot of differentially expressed genes between CA1 and CA3 clusters. b. Heatmap of the expression level of all genes between the two groups. c. CHOIR-based testing between each reported cluster and its nearest neighbor. Valid means the cluster is significantly different from its neighboring cluster and should be split from each other, while invalid means two clusters are not significantly different and should be merged.\n(2) Additionally, no baseline methods are included for comparison on this dataset, making it difficult to evaluate the distinct advantages or effectiveness of RareQ in this context.\nResponse: We analyzed the dataset using several baseline methods—scCAD, CellSIUS, FiRE, and GapClust—for comparison. Neither scCAD nor CellSIUS detected any rare cell clusters. Although FiRE and GapClust identified some rare clusters, their spatial distributions were highly dispersed and lacked meaningful anatomical structure (Fig. R15a, b). These findings suggest that existing methods struggle to reliably identify rare cell populations and assign them to their correct biological locations. In contrast, RareQ accurately identified and localized rare cells within the rostral migratory stream, the hippocampus, and the Purkinje layer of the cerebellum. This demonstrates RareQ’s distinct advantage in both detecting and spatially resolving rare cell populations within their proper anatomical context.\nThese comparative results have been incorporated into the revised manuscript.\nFig. R15. Rare cell types predicted by FiRE (a) and GapClust (b) in EEL FISH mouse brain data\n8. The discovery of rare cell types in Figure 7d appears to be driven largely by over-clustering. In the Xenium dataset, 162,000 cells were divided into 199 clusters—most containing fewer than 1% of cells. This suggests that the identification of rare populations stems more from arbitrary clustering granularity than from any unique strength of the method. Standard clustering tools could yield similar results, and the approach is highly sensitive to minor cellular heterogeneity, casting serious doubt on the biological validity of these rare cell types.\nResponse: To ensure that RareQ’s identification of rare populations was not an artifact of over-clustering, we applied CHOIR (Cluster Hierarchy Optimization by Iterative Random forests) (https://doi.org/10.1038/s41588-025-02148-8), a recently published statistical framework specifically designed to guard against both over- and under-clustering in single-cell data. Using CHOIR’s compareClusters function, we evaluated the 199 clusters identified by RareQ in the Xenium dataset. CHOIR confirmed that 189 of these clusters represented statistically distinct populations (Fig. R16a), indicating that the vast majority of RareQ-defined clusters are biologically meaningful rather than spurious subdivisions.\nImportantly, RareQ is not simply a high-resolution clustering approach; it leverages local neighborhood connectivity (Q) to distinguish genuine rare populations from arbitrary splits. To further validate this, we compared RareQ with Seurat’s Louvain clustering across a wide range of resolutions (0.8, 1, 2, 5, and 10). At default and moderate resolutions (0.8–2), Seurat failed to detect the CA2 region (Fig. R16b). At higher resolutions (5 and 10), Seurat recovered CA2 but was unable to resolve the CA3_1 and CA3_2 subpopulations, which RareQ successfully distinguished.\nTogether, these results—combining CHOIR’s statistical validation with comparative benchmarking—demonstrate that RareQ can identify rare populations with high specificity while maintaining biological plausibility.\nFig. R16 a, CHOIR-based testing between each cluster and its nearest neighbor predicted by RareQ in Xenium mouse brain data. Valid means the cluster is significantly different from its neighboring cluster and should be split from each other, while invalid means two clusters are not significantly different and should be merged. b, Spatial maps of Xenium mouse brain data, where every dot is a single cell colored by clusters predicted by Seurat at specified resolutions and a magnified view showing the hippocampus labeled by identified clusters.\n9. The description of the methodology on page 13, lines 468–477, is notably unclear and problematic.\n(1)The parameter k mentioned here should be clearly distinguished from the k value defined in the Data Preprocessing section to avoid confusion.\nResponse: We appreciate the reviewer’s observation. In the Data Preprocessing section, the k value refers to the parameter used for constructing the kNN graph, whereas in RareQ, the k value denotes the parameter used in calculating the Q value. To avoid confusion, we have renamed the former as k.param and added further clarification to clearly distinguish between these two definitions.\n(2)According to the manuscript’s logic, since all cells in the set S are connected as nearest neighbors, this would imply that every cell is connected by edges, resulting in all calculated Q values being equal to 1.\nResponse: We thank the reviewer for this comment. It's important to clarify that Q will equal 1 only when a cell and its k-nearest neighbor network form a perfect cliquishness, meaning every member is connected exclusively to every other member in the set S (Fig. R17a). However, when there exist edges from cells in S to external cells (i.e., outward connection (OC) > 0) (Fig. R17b, c), the Q value would be less than 1, as the denominator |IC| + |OC| is always greater than the numerator |IC|.\nFig. R17 Schematic plot of three different connectivity scenarios of query cell’s neighborhood with different Q values.\n10. Accurate identification of rare cells should be based on a correct overall clustering performance. However, the manuscript lacks evaluation metrics assessing the quality of the global clustering results.\nResponse: We thank the reviewer for raising this important point. In our initial submission, we primarily evaluated performance using metrics widely applied in published rare cell detection studies, such as F score. While we did assess global clustering quality on the PBMC simulation scRNA-seq datasets (Fig. R18a), we acknowledge that this evaluation was not extended to all datasets.\nTo address this, we have now systematically evaluated the overall clustering performance using the Normalized Mutual Information (NMI) metric across a broader range of datasets. On one hand, we have benchmarked RareQ against the clustering algorithms (CellSIUS, RaceID, GiniClust2, EDGE) tailored for scRNA-seq dataset in the four PBMC scRNA-seq datasets and 20 real scRNA-seq datasets (Fig. R18b, c). On the other hand, we compared RareQ with multi-omics approach MarsGT in the 150 paired PBMC simulation datasets, 4 PBMC multiome datasets, and 10 paired scRNA-seq and scATAC-seq datasets (Fig. R19a, b and Fig. R11d). As shown in the new results, RareQ consistently achieved the highest NMI scores compared with other algorithms, demonstrating both accurate rare cell identification and robust global clustering performance. These results have been incorporated into the revised manuscript (Supplementary Figure 22, 46).\nFig. R18. Benchmarking RareQ against existing methods on global clustering performance in simulated and real scRNA-seq datasets. a. NMI values evaluated across 150 PBMC simulated scRNA-seq datasets by specified methods. b. Boxplots showing NMI values evaluated across four PBMC scRNA-seq datasets by specified methods. c. Boxplots showing NMI values evaluated across 20 real scRNA-seq datasets by specified methods.\nFig. R19. Benchmarking RareQ against MarsGT on global clustering performance in simulated and real paired scRNA-seq and scATAC-seq datasets. a. NMI values evaluated across 150 paired PBMC simulated datasets by specified methods. b. Boxplots showing NMI values evaluated across four paired PBMC datasets by specified methods.\n11. Given the discussion on memory usage and runtime, it is necessary for the authors to specify the exact GPU and CPU models used in their experiments.\nResponse: Thank you for this helpful suggestion. We have added detailed hardware specifications to the revised manuscript. In most experiments, we used a computing platform with 64 CPUs (Intel® Xeon® Gold 6226R, 2.90 GHz), 750 GB of memory, and 2 GPUs (NVIDIA RTX A6000, 48 GB memory each). For benchmarking MarsGT, we employed a platform with 128 CPUs (AMD EPYC 7543, 32-Core Processor).\n12. Since the study only uses spatial transcriptomics data for testing, the umbrella term “spatial omics” is misleading and should not be used. The authors should add at least one additional spatial modality, for example, spatial ATAC-seq, to substantiate multi-omic applicability.\nResponse: We thank the reviewer for raising this point. To address the concern, we extended our analysis to an in situ gene–protein expression dataset derived from a human renal cell carcinoma FFPE sample by the Xenium platform.\nUsing RareQ, we identified 59 distinct cell clusters forming complex spatial architectures (Fig. R6a). These clusters were clearly separated in UMAP space, reflecting unique transcriptional and proteomic signatures, including major populations such as B cells, T cells, myeloid cells, tumor cells, fibroblasts, and endothelial cells (Fig. R6b, c).\nRareQ further uncovered several biologically meaningful features:\n­ Tertiary lymphoid structures (TLS): In the peritumoral region, RareQ resolved highly organized aggregates of B cells consistent with TLS (Fig. R7a).\n­ Myeloid heterogeneity: RareQ distinguished diverse myeloid subpopulations, including mast cells, conventional dendritic cells (cDCs), plasmacytoid dendritic cells (pDCs), mature cDCs, and monocytes, each defined by distinct gene- and protein-level markers (Fig. R7b).\n­ Rare proliferative subclusters: RareQ detected rare Ki-67⁺ proliferative clusters (clusters 1, 5, 7, 12, and 14) across tumor, fibroblast, CD4⁺ T-cell, and plasma cell lineages (Fig. R7c, d).\nIn contrast, methods such as CellSIUS, FiRE, and GapClust—when applied to gene expression data alone—failed to resolve subclusters with comparable resolution (Fig. R8).\nTogether, these results demonstrate RareQ’s capacity to analyze spatial multi-omic datasets with high precision.\n13. In line 468, the parameter k is set directly to 6. Since k is a critical parameter, its impact on the experimental results should be explicitly discussed, with a dedicated analysis of how varying k affects performance.\nResponse: We thank the reviewer for this valuable suggestion. To address this, we conducted a comprehensive sensitivity analysis by systematically varying k from 5 to 30 across 150 simulated and 20 real scRNA-seq datasets.\nOur analysis showed that RareQ's performance is highly robust to changes in k. In the simulated PBMC datasets, both the F scores (measuring rare cell detection accuracy) and NMI values (measuring overall clustering quality) remained stable with only minimal fluctuations (Fig. R20a, b). We observed the same trend in the 20 real datasets, where the F , precision, recall, and NMI values consistently stayed high across the entire range of k values (Fig. R20c, d).\nThese results demonstrate that RareQ is not overly sensitive to the choice of k, and its performance remains stable across a wide range of values. We have included these findings in the revised manuscript.\nFig. R20. Evaluation of RareQ’s sensitivity to the choice of k value. a. Mean F scores (top) and detailed F scores (bottom) across different simulation datasets with k values ranging from 5 to 30. b. NMI values evaluated across the same simulation datasets and k ranges as in panel (a). c&d. Boxplots showing the distribution of F score, precision, recall, and NMI at varying k values across 20 real scRNA-seq datasets.\n14. The absence of example code in the provided code repository hinders reproducibility and makes it difficult for users to effectively apply the method.\nResponse: We have updated the code repository (www.github.com/xiaolab-xjtu/RareQ) to include example code, sample datasets, and a step-by-step tutorial."}, {"round": 2, "reviewer_comments": {"reviewer_1": "Thank you for the revision. Despite the additional datasets and sensitivity analyses, the revision does not resolve several fundamental concerns raised in my original report, and the manuscript still falls short of the conceptual, methodological, and biological standards for Nature Communications.\n1. The “Q value” corresponds to the classical clustering coefficient and reframes the work as an application of this index to rare cell detection. However, beyond adopting this metric, the manuscript still does not clearly articulate what is methodologically new relative to standard graph based clustering/label propagation pipelines (e.g., how RareQ’s propagation/merging is fundamentally different from existing deterministic label propagation plus post hoc merging). Please provide the algorithmic innovations with ablations isolating the contribution of each step and benchmark these ablations directly.\n2. The definition is “low abundance and biological significance evidenced by distinct molecular profiles”. The supporting evidence largely relies on DEGs, GSEA, and signature scores, which are circular with respect to the same transcriptomic features used for discovery. No orthogonal validation (imaging, flow, perturbation, or independent pathological assessment) is provided.\n3. The added lymphoma, AD snRNA seq, and RCC spatial multi omics datasets are welcome. However, the analyses do not demonstrate control for donor level effects, batch, or case/control imbalance, particularly for the AD microglia where cluster 15 is “enriched in AD”. What is the statistical framework used to test group specific enrichment after accounting for donor as a random effect, and potential covariates; otherwise the disease association claims are difficult to interpret.\n4. One of my important concerns was that RareQ’s annotations on public spatial datasets did not match the original study reports. I do not see a systematic, side by side reconciliation with the original labels or a documented re annotation protocol.\n5. CHOIR is applied asymmetrically, with no equivalent CHOIR validation for high resolution Seurat baselines under comparable features and parameters. Moreover, because CHOIR builds RF classifiers on the same expression space used to define clusters, there is a risk of self validation.\n6. Were MarsGT hyperparameters tuned comparably? Were identical kNN graphs/latent spaces supplied where applicable? Were input features harmonized (e.g., WNN vs. integrated embeddings)?\n7. In response to my request to avoid the umbrella term “spatial omics” if only RNA was used, you added one Xenium RNA+ADT dataset. This is helpful but still limited to a single sample and modality pair. Either temper the general claim (e.g., “spatial transcriptomics and one transcript protein dataset”) or broaden the evidence (additional tissues/samples or orthogonal spatial modalities).\n8. The tie break rule (smallest index) is deterministic but the biological stability of final clusters under index permutations is only assessed via aggregate metrics (F1/NMI).\n9. In challenging spatial settings, RareQ may not add biological insight beyond recapitulating known structures.\n10. To facilitate verification, please (i) include end to end scripts that reproduce a subset of the new disease analyses (e.g., AD microglia and RCC), (ii) pin package versions and random seeds, and (iii) provide exact commands for all figures cited in the response (R series).\nAfter carefully reviewing the point‑by‑point response and new analyses, I find that the core concerns remain unresolved. In its current form, the manuscript is not suitable for publication in Nature Communications.", "reviewer_2": "The authors have responded comprehensively to my comments.\nThe additional analyses would justify the choice of hyperparameters and should be included in the final manuscript.", "reviewer_3": "My comments have been properly addressed. I recommend the publication of the manuscript. (Zixuan Cang)"}, "author_response": "We thank the reviewer for the rigorous standards and constructive critique, which was essential for enhancing the quality of our manuscript. We have made extensive revisions to improve the methodological clarity of RareQ and to more effectively illustrate its capacity to yield novel scientific insights. These include ablation studies to elucidate the algorithmic innovations, integration of additional biological validation datasets, and implementation of reproducible analysis scripts to ensure transparency and replicability. We are confident that these changes now align with the rigorous criteria for publication in Nature Communications.\n1. The “Q value” corresponds to the classical clustering coefficient and reframes the work as an application of this index to rare cell detection. However, beyond adopting this metric, the manuscript still does not clearly articulate what is methodologically new relative to standard graph based clustering/label propagation pipelines (e.g., how RareQ’s propagation/merging is fundamentally different from existing deterministic label propagation plus post hoc merging). Please provide the algorithmic innovations with ablations isolating the contribution of each step and benchmark these ablations directly.\nResponse:\nWe thank the reviewer for raising these concerns regarding the methodological novelty of RareQ and for providing the opportunity for further clarification. While the Q value is indeed inspired by the classical clustering coefficient in small-world network theory, RareQ transforms this concept into a dynamic, cell-specific topological measure that quantifies the interconnectivity within each cell’s neighborhood network. It is also extended to the cluster level (denoted Qc) to guide post-hoc merging for optimizing clustering resolution and preventing over-clustering.\nIt is noteworthy that several established rare-cell identification algorithms—such as GiniClust, RaceID, and FiRE—also adapt well-established concepts to the single-cell context, such as Gini impurity, negative binomial-based outlier modeling, and sketching with locality-sensitive hashing. The innovation of such methods—including RareQ—lies not in inventing a new mathematical index, but in its algorithmic operationalization to solve specific tasks. RareQ’s directed label propagation based on Q value is specifically tailed for identifying small and topologically self-organized subpopulations, which differs from conventional approaches based on uniform diffusion or modularity optimization (e.g., Louvain, Leiden). Below, we summarize RareQ’s key algorithmic innovations and ablation studies for benchmarking the contribution of its key components.\n(i) Quantitative characterization of rare-cell neighborhood topology\nRareQ distinguishes itself from local density estimation methods by characterizing the topological structure of the cell neighborhood graph; specifically, it quantifies the degree of interconnectivity among nodes within each cell’s k-nearest-neighbor graph. This enables identification of tightly connected cell cliques. Cells within smaller clusters exhibit higher Q values compared to those in larger populations, highlighting the effectiveness of the Q value in characterizing their relative topological isolation, as demonstrated in Figure 1b and Figure 5b of the manuscript.\n(ii) Q-guided directed label propagation\nRareQ employs a Q-guided label-propagation strategy that enables the detection of rare cell clusters. This includes not only topologically isolated groups but also, crucially, those embedded within or adjacent to major populations—a challenge for conventional graph-based or density-based frameworks where such rare cells are often absorbed into their abundant neighbors. RareQ overcomes this limitation by introducing a cell-specific Q metric that quantifies local topological coherence. During label propagation, cells with lower Q values adopt the label of the highest-Q neighbor within their local neighborhood. In this way, high-Q rare-cell cliques function as stable, topology-driven anchors that guide label diffusion outward, establishing a directional propagation process. This mechanism effectively captures both types of challenging rare populations. To assess the importance of this design, we performed ablation studies in which the Q-guided propagation was replaced by uniform diffusion (equivalent to assigning a constant Q = 0). This substitution resulted in decreased F1 scores, precision, and recall for several rare clusters (Fig. R1a–e), underscoring the essential role of Q-guided propagation in accurately identifying rare cell populations.\nWe speculated that this decline specifically affects clusters that have weak separation from dominant populations, as those completely isolated within the KNN graph would be easily detected by conventional methods. To test our hypothesis, we created detection difficulty through two distinct approaches: (1) we subsampled rare cells from major clusters and perturbed the expression of 10–100 genes by three-fold up regulation (Fig. R1f). (2) we adapted the real data shown in Supplementary Figure 23 and incorporated varying numbers of differentially expressed genes into the real rare clusters (Fig. R1g). The results showed that RareQ consistently outperformed the ablated version, confirming that Q-guided propagation enhances identification of challenging populations. Moreover, in practical applications such as the AD snRNA-seq dataset, several clusters became undetectable when the Q-values were removed (Fig. R2a, b).\nFig. R1. Ablation analyses of Q-value demonstrate its contribution to rare-cell identification across simulation and real datasets\na. Cumulative curves comparing F1 scores before and after Q-value ablation in simulated datasets 1 and 2.\nb. Same as (a), but based on simulated dataset 3.\nc. Same as (a), but for NMI comparison in three simulated datasets.\nd. Comparison of F1 score, precision, and recall before and after Q-value ablation in real datasets.\ne. Same as (d), but for NMI comparison in the real datasets.\nf. Ablation analysis of the impact of Q-value under simulated conditions with varying degrees of gene-expression perturbations.\ng. Same as (f), but examining the impact of varying the number of differentially expressed genes in the simulation.\nFig. R2. Ablation analysis of the Q-value demonstrates its contribution to rare-cell identification in the Alzheimer’s disease dataset\na. UMAP plots showing the clusters before and after Q-value ablation in Alzheimer’s disease (AD) snRNA-seq dataset.\nb. Confusion matrix comparing the identified cluster before and after Q ablation.\n(iii) ΔQc-based topology-aware recursive merging\nFollowing propagation, RareQ employs a recursive merging strategy driven by the cluster-level Q metric (Qc). Merges are accepted only if the resulting cluster exhibits a higher Qc value, reflecting increased internal connectivity. This provides a biologically interpretable, topology-based criterion that prevents over-clustering while preserving stable rare-cell clusters. Ablation analysis omitting this step led to a slightly increased F1 score but significantly decreased NMI, confirming that ΔQc-guided merging significantly improves overall global clustering accuracy, despite potential underclustering in a few clusters (Fig. R3).\nFig. R3. Ablation analysis of the merging component demonstrates its contribution to guarantee the clustering quality\na. Cumulative curves comparing F scores before and after merging step ablation in 1 simulated datasets 1 and 2.\nb. Same as (a), but based on simulated dataset 3.\nc. Cumulative curves comparing F scores before and after merging ablation in 20 real scRNA-seq datasets.\nd. Same as (c), but comparing NMI.\nIn summary, RareQ's core methodological innovation is the extension of a static small-world descriptor into a topology-aware propagation and merging strategy. This approach significantly enhances the detection of hard-to-detect rare cells while maintaining the integrity of major clusters. The corresponding ablation analyses are presented in supplementary figure 65-66, and all scripts are openly accessible on GitHub (https://github.com/xiaolab-xjtu/RareQ-reproducible-scripts) to ensure reproducibility.\n2. The definition is “low abundance and biological significance evidenced by distinct molecular profiles”. The supporting evidence largely relies on DEGs, GSEA, and signature scores, which are circular with respect to the same transcriptomic features used for discovery. No orthogonal validation (imaging, flow, perturbation, or independent pathological assessment) is provided.\nResponse:\nWe appreciate the reviewer revisiting this critical concern and the opportunity for further clarification. We acknowledge the reviewer’s concern that differential gene expression (DEG) analysis, gene set enrichment analysis (GSEA), and signature scoring might appear circular, as their purpose was to evaluate whether the identified rare populations exhibit coherent and biologically interpretable molecular programs.\nWe would like to reiterate that, for the Xenium mouse brain spatial transcriptomics dataset, we have already provided extensive orthogonal evidence in the previous revision to substantiate the rare CA3 subpopulations identified by RareQ. The validation data include in situ hybridization, RNAscope, riboprobe-based assays, and an independent scRNA-seq dataset. These data consistently confirmed the specificity and existence of the CA3 subpopulations specifically detected by RareQ. In this revision, we further validated our findings using two independent mouse brain spatial transcriptomics datasets from 10x Genomics Datasets (https://www.10xgenomics.com/datasets) (Fig. R4).\nFig. R4. Independent spatial transcriptomics datasets validate the spatial localization of CA3 subpopulations and their marker gene’s expression.\na. Expression patterns of marker genes for CA2, CA3, and the two CA3 subtypes in the mouse brain spatial transcriptomics dataset.\nb. Same as (a), showing the expression patterns in another mouse brain dataset.\nMarkers of CA3_1: Strip2; CA3_2: Rspo2, Plcxd2, Nwd2.\nTo achieve the most rigorous validation possible, we have now incorporated several additional validation datasets, specifically:\n(1) For the airway epithelial dataset, we leveraged three independent human spatial transcriptomics datasets (Murthy P. et al., 2022, PMID: 35355018; Madissoon E. et al., 2023, PMID: 36543915) and confirmed the presence of the RareQ-discovered Cycling Basal (C3), Ciliated (C4), and Secretory cell (C10) populations in specific regions along the airway lumen (Fig. R5), as defined by their RareQ-identified markers (Fig. 5d, e). These spatial datasets provide both molecular and contextual information that surpasses the limitations of bulk analysis (like flow cytometry) or purely morphological imaging approaches, thereby enabling significantly more precise validation.\nFig. R5. Validation of RareQ's unique rare airway epithelial populations using independent spatial transcriptomics.\nSpatial distribution of spots highly expressing marker genes of C3/C6 (a), C4 (b), and C10 (c).\nMarkers of C3: TOP2A– MKI67+ H2AFZ+ KRT5+ TP63+ NGFR+; C6: TOP2A+ MKI67+ H2AFZ+ KRT5+ TP63+ NGFR+; C4: CCNO+ FOXJ1+ LRRC23+ MEIG1+ TPPP3+; C10: SCGB1A1+ FOXJ1+ LRRC23+ TUBA1A+ ISG15– IFIT3–.\n(2) For the Alzheimer’s disease (AD) dataset, we validated the significant enrichment of the C15 microglial population signature in AD patients compared with healthy controls using an independent bulk RNA-seq cohort (Fig. R6).\nFig. R6. Independent bulk RNA-seq confirms the significant enrichment of C15 microglial populations in Alzheimer’s disease (AD) patients versus control individuals.\n(3) In the B cell lymphoma datasets, the simultaneous scRNA-seq and scATAC-seq inherently provide cross-modality data validation. Furthermore, rare cell types, such as the PD-L1+ dendritic cells specifically identified by RareQ, have been independently reported in other published studies (Zilionis R et al., 2019, PMID: 30979687; Zhang Q et al., 2019, PMID: 31675496; Maier, B. et al., 2020, PMID: 32269339).\nCollectively, these multi-level validations provide independent, reproducible, and biologically grounded support for the rare populations detected by RareQ. To clarify, performing additional orthogonal experimental validations for every rare cell cluster would require specialized biological systems and substantial resources that fall beyond the scope of this computational methodology study. Nevertheless, our current findings establish a solid computational and biological foundation for future dedicated experimental investigations.\n3. The added lymphoma, AD snRNA seq, and RCC spatial multi omics datasets are welcome. However, the analyses do not demonstrate control for donor level effects, batch, or case/control imbalance, particularly for the AD microglia where cluster 15 is “enriched in AD”. What is the statistical framework used to test group specific enrichment after accounting for donor as a random effect, and potential covariates; otherwise the disease association claims are difficult to interpret.\nResponse: We thank the reviewer for this insightful comment. To account for donor-level effects, we performed differential cell abundance analysis using a negative binomial generalized linear model (NB-GLM) with quasi-likelihood (QL) estimation, consistent with the approach implemented in the Milo algorithm (Dann E et al., 2022, PMID: 34594043). In this model, donor identity was included as a random effect within a quasi-likelihood mixed model to account for intra-donor correlations and inter-donor variability. Group-specific enrichment (e.g., Alzheimer’s disease versus control) was then assessed using a QL F-test with specified contrasts for each cluster.\nAs shown in the Table R1 and Supplementary Table 6, among the three microglial subtypes (clusters 15, 16, and 21), cluster 15 exhibited a statistically significant enrichment of AD-derived cells after adjusting for donor effects, consistent with previously reported AD-associated microglial phenotypes. This statistical framework has now been explicitly detailed in the revised Methods section.\nTable R1. Differential abundance of cell subtypes in AD vs. healthy controls\nCluster Log2 FC LogCPM F-statistic p value\n1 -1.5403 10.1705 1.0636 0.3040\n2 3.6042 10.0398 4.0582 0.0457\n3 0.8819 10.6279 0.3307 0.5660\n4 2.0888 11.4919 1.1927 0.2765\n5 5.1107 11.5020 4.3424 0.0388\n6 -0.0828 11.9391 0.0028 0.9573\n7 -2.9078 12.1873 3.2537 0.0732\n8 -1.5442 12.4684 1.0134 0.3156\n9 0.0329 11.6586 0.0003 0.9848\n10 -0.0552 12.4467 0.0014 0.9701\n11 -3.4253 12.5718 4.1951 0.0422\n12 1.0213 12.0423 0.4811 0.4889\n13 -0.2559 13.0149 0.0393 0.8429\n14 0.7780 12.6065 0.3025 0.5830\n15 3.5935 13.1931 5.4319 0.0211\n16 0.0236 13.1318 0.0003 0.9852\n17 5.1205 13.4050 8.8075 0.0034\n18 -2.2309 14.7938 2.1299 0.1465\n19 -4.5116 14.8527 6.5940 0.0112\n20 4.0969 14.7445 5.3698 0.0218\n21 0.2677 15.7650 0.0533 0.8177\n22 0.7794 14.9526 0.3262 0.5687\n23 -0.5291 16.2058 0.2086 0.6484\n24 -0.0865 17.4168 0.0055 0.9413\n25 0.3220 19.3085 0.3615 0.5485\n4. One of my important concerns was that RareQ’s annotations on public spatial datasets did not match the original study reports. I do not see a systematic, side by side reconciliation with the original labels or a documented re annotation protocol.\nResponse: We thank the reviewer for raising this important point. Among the three public spatial datasets analyzed in our study, only the EEL FISH dataset included cell-type annotations in the original publication; the remaining datasets were released by the 10x Genomics platform without accompanying cell-type labels. As explained in our previous revision, for the EEL FISH dataset, we were unable to reproduce the original annotations even when re-clustering with more sensitive parameters following the given methods. Given this limitation, we have moved the results of this dataset to the Supplementary Figure 54.\nIn this revision, we conducted a systematic, side-by-side reconciliation of RareQ-derived annotations with the original cluster labels for the EEL FISH dataset, which is presented in the newly added Supplementary Figure 54d. This analysis revealed that the original annotations contained substantially more clusters than those identified by RareQ (Fig. R7a). We speculate that, if spatial context is not taken into account, the original annotations based solely on expression profiles of 440 genes may have been over-clustered, as demonstrated in our response to Comment 5. We also re-annotated this dataset following the procedures described in the original publication, using a resolution that produced a comparable number of clusters, and obtained similar results (Fig. R7b). To enhance transparency and reproducibility, we have provided a comprehensive re-annotation protocol on Github (https://github.com/xiaolab-xjtu/RareQ-reproducible-scripts).\nFig. R7. Confusion matrices comparing RareQ clusters with reference annotations.\na. RareQ vs. original annotations.\nb. RareQ vs. Seurat re-annotations.\n5. CHOIR is applied asymmetrically, with no equivalent CHOIR validation for high resolution Seurat baselines under comparable features and parameters. Moreover, because CHOIR builds RF classifiers on the same expression space used to define clusters, there is a risk of self validation.\nResponse: We have now conducted CHOIR analyses under comparable settings for both RareQ and the high-resolution Seurat baselines. In the Xenium spatial transcriptomics dataset, CHOIR indicated that most clusters identified by both RareQ and Seurat showed low probabilities of over-clustering (Fig. R8). However, for the EEL FISH spatial dataset, CHOIR evaluation suggested that the majority of clusters from the original annotations (116 of 188) and from Seurat re-clustering (73 of 188) were flagged as over-clustered (Fig. R9), only a minor subset of RareQ-derived clusters exhibited this issue (Supplementary figure 57a).\nFig. R8. Assessment of over-clustering using CHOIR for RareQ-derived clusters (a) and Seurat-derived clusters (b) on Xenium spatial transcriptomic data.\nFig. R9. Assessment of over-clustering using CHOIR for original clusters (a) and Seurat reannotated clusters (b) on EEL Fish spatial transcriptomic data.\nWe appreciate the reviewer's concern regarding potential self-validation within the CHOIR analysis, as the classifier operates on the same expression space used for clustering. It's important to note that RareQ inherently mitigates over-clustering through its dedicated post hoc merging module. Furthermore, we used the Silhouette Coefficient—an effective internal metric for assessing cluster count appropriateness and overall quality—to independently evaluate RareQ’s clustering quality against Seurat (Fig. R10a). This analysis showed that RareQ's clusters generally achieved a higher quality score than those generated by Seurat. For the separate EEL FISH dataset, we further demonstrated that clusters derived from both the original annotation and Seurat re-clustering were generally of low quality (Fig. R10b).\nFig. R10. Assessment of cluster quality using silhouette coefficients\na. Assessment of RareQ-derived clusters and Seurat-derived clusters on Xenium spatial transcriptomic data.\nb. Assessment of original clusters and Seurat reannotated clusters on EEL Fish spatial transcriptomic data.\nFurthermore, the RareQ's unique rare clusters showed distinct biological identities, evidenced by differential gene expression when compared to their most similar major clusters (Fig. R11). Finally, as detailed in our response to Comment 2, we conducted orthogonal validation of the RareQ-identified clusters. Across all analyses performed, RareQ consistently demonstrated that the resulting clusters are biologically meaningful and not artifacts of over-clustering.\nFig. R11. Heatmap illustrating the molecular differences between RareQ-specific rare cell clusters and their similar major cell types.\n6. Were MarsGT hyperparameters tuned comparably? Were identical kNN graphs/latent spaces supplied where applicable? Were input features harmonized (e.g., WNN vs. integrated embeddings)?\nResponse: Following the original paper's methodology, we performed a grid search using four out of the ten multiome datasets to optimize MarsGT's hyperparameter combination, including weighted decay, learning rate, and label smoothing. Our results indicated that MarsGT achieves its best performance with the settings: weighted decay (wd) = 0.3, learning rate (lr) = 0.0005, and label smoothing (labsm) = 0.3 (Fig. R12a). This set was subsequently adopted as the optimized parameter combination for all ten multiome datasets.\nWe further input the WNN-integrated features into MarsGT (denoted as MarsGT_WNN), matching the input used for RareQ. While MarsGT_WNN showed substantial improvements in rare-cell detection accuracy (F score, precision, and recall) and overall clustering performance (NMI) compared to the default MarsGT input, its rare-cell detection accuracy remained slightly lower than that of RareQ_WNN (Fig. R12b). This new comparative analysis has been incorporated into the revised manuscript, and the results are discussed. Note that due to the long runtime of MarsGT, it is difficult to complete testing on large-scale simulated datasets. Therefore, we did not include MarsGT_WNN results on the simulation datasets.\nFig. R12. Hyperparameter tuning of MarsGT and performance comparison under its optimal parameter settings.\na. F score comparison across different combinations of weight decay, learning rate, and label-smoothing parameters.\nb. Performance comparison among MarsGT results obtained using the optimal parameter set and WNN-integrated embeddings as input (MarsGT_WNN) and those generated by RareQ.\n7. In response to my request to avoid the umbrella term “spatial omics” if only RNA was used, you added one Xenium RNA+ADT dataset. This is helpful but still limited to a single sample and modality pair. Either temper the general claim (e.g., “spatial transcriptomics and one transcript protein dataset”) or broaden the evidence (additional tissues/samples or orthogonal spatial modalities).\nResponse: We thank the reviewer for this thoughtful comment. In the revised manuscript, we now describe the application of RareQ to spatial transcriptomics datasets and one transcript–protein (RNA+ADT) dataset, replacing the broader term “spatial omics.” This revision provides a more accurate and precise description of the modalities analyzed while maintaining clarity.\n8. The tie break rule (smallest index) is deterministic but the biological stability of final clusters under index permutations is only assessed via aggregate metrics (F1/NMI).\nResponse: In response, we further evaluated clustering stability using the Jaccard index, a non-aggregate metric sensitive to local cell label differences. This analysis assessed the consistency of both cell memberships and cluster-specific gene signatures across runs. The results showed that random index selection produced outcomes nearly identical to those of the deterministic rule, with only minor variance observed in a few simulations (Fig. R13). To enhance flexibility, we have also implemented an optional non-deterministic tie-break mode, allowing users to perform multiple random runs and derive consensus clusters for obtaining higher-precision cluster assignments.\nFig. R13. Evaluating the impact of tie-break rules on RareQ stability\na. Effect of index permutation on the consistency of cell membership assignments in three simulation datasets.\nb. Effect of index permutation on the stability of the top 50 highly specific marker genes per cluster in the tree simulation datasets.\nc. similar to (a) but test on the 20 real datasets.\nd. similar to (b) but test on the 20 real datasets.\n9. In challenging spatial settings, RareQ may not add biological insight beyond recapitulating known structures.\nResponse: We thank the reviewer's concern. Spatial transcriptomics data are conceptually similar to scRNA-seq data but incorporate spatial information. In our study, RareQ was applied using expression data only, without incorporating spatial priors, making its accuracy in spatial transcriptomic analysis analogous to that in scRNA-seq. Therefore, we anticipate that RareQ can provide comparable biological insights in spatial contexts, as demonstrated in the Xenium mouse brain data, where RareQ uncovered CA3 subpopulations within known anatomical regions that competing methods failed to detect.\nIn this revision, we extend these findings by uncovering additional biological insights:\nIn the EEL FISH mouse brain dataset, RareQ identified a rare choroid plexus epithelial subpopulation (C1) characterized by high expression of Foxj1, Otx2, and Lratd2, spatially nested within the dominant epithelial cluster (C11) (Fig. R14a-c). The spatial localization of the C1 population, along with its high expression of cilium-associated genes such as Foxj1, suggests that this group likely corresponds to ciliated epithelial cells of the choroid plexus. Previous studies have shown that cilia play a critical role in regulating cerebrospinal fluid flow and maintaining the ventricular microenvironment. Therefore, the rare ciliated-like choroid plexus epithelial population identified in our study provides a higher-resolution cellular reference for elucidating the fine structural organization and functional architecture of the ventricular system. Consistently, in another independent spatial dataset from an early developmental stage of the mouse brain, we observed a corresponding epithelial population marked by Foxj1 and Otx2 positioned adjacent to the main choroid plexus clusters defined by Ttr (Fig. R14d, e). We further validated these findings using an independent, high-resolution choroid plexus scRNA-seq dataset, where epithelial cells could be precisely mapped to the corresponding clusters in the EEL FISH spatial data (Fig. R14f, g). Of the 16 epithelial populations detected by RareQ, clusters like C3, C4, and C5 showed elevated expression of Foxj1, Otx2, and Lratd2, mirroring the characteristics of the spatially-defined C1 cluster (Fig. R14h). Together, these observations indicate that the RareQ identified rare cells in spatial context represent a biologically meaningful, specialized functional state, rather than technical noise.\nOverall, RareQ not only recapitulates known spatial structures but also detects subtle and biologically meaningful rare-cell populations.\nFig. R14. Identification and characterization of choroid plexus cell subtypes.\na. Original annotated spatial map of the choroid plexus population.\nb. Spatial distribution of the C1 and C11 choroid plexus populations identified by RareQ.\nc. Bubble plot showing marker genes of the C1 population.\nd. Spatial distribution of diverse cell types of a mouse head sagittal section.\ne. Expression patterns and spatial localization of the choroid plexus marker genes.\nf. RareQ-identified choroid plexus cell populations in mouse brain snRNA-seq data.\ng. Cross-modal correlation of epithelial cells (snRNA-seq) with all cells in the EEL FISH spatial dataset.\nh. Bubble plot illustrating marker genes for the c1 choroid plexus population within the snRNA-seq cell clusters.\n10. To facilitate verification, please (i) include end to end scripts that reproduce a subset of the new disease analyses (e.g., AD microglia and RCC), (ii) pin package versions and random seeds, and (iii) provide exact commands for all figures cited in the response (R series).\nResponse: We have implemented four tutorials and supplied the requested end-to-end scripts to guarantee complete reproducibility at GitHub repository (https://github.com/xiaolab-xjtu/RareQ/tree/main/tutorials; https://github.com/xiaolab-xjtu/RareQ-reproducible-scripts). Specifically, this repository includes:\n1) Complete scripts for all new analyses, including the Alzheimer’s disease (AD) microglia and renal cell carcinoma (RCC) case studies.\n2) Pinned package versions and dependencies, with fixed random seeds utilized across all analyses.\n3) A detailed guide containing the exact R commands required to reproduce every figure cited in our responses.\nAfter carefully reviewing the point-by-point response and new analyses, I find that the core concerns remain unresolved. In its current form, the manuscript is not suitable for publication in Nature Communications.\nResponse: We appreciate the reviewer’s continued engagement and careful evaluation. In this revision, we have thoroughly addressed all concerns through extensive new analyses, expanded validation datasets, and detailed methodological clarifications. We respectfully hope the reviewer will reconsider the revised manuscript in light of the new evidence and improvements provided."}, {"round": 3, "reviewer_comments": {"reviewer_1": "This revised manuscript introduces RareQ, a graph-based clustering framework designed for rare-population discovery through Q-guided label propagation and Qc-guided recursive merging. While the authors have addressed some previous comments, several fundamental methodological concerns remain unresolved. Most notably, the justification for the method's performance and hyperparameter selection relies heavily on supervised metrics, which contradicts the \"unsupervised discovery\" premise. Furthermore, the distinction between RareQ and existing graph-clustering pipelines remains insufficiently demonstrated.\n1. If RareQ is intended for unsupervised discovery, the reliance on label-dependent metrics (e.g., F1-score, NMI, Precision/Recall) to justify hyperparameter choices and \"high-accuracy regions\" is problematic. These metrics are unavailable in real-world discovery scenarios. Specifically, the suggestion in Supplementary Note 2 that users \"modify parameters accordingly\" to improve accuracy is impractical without ground-truth labels. The authors must provide a fully label-free parameter selection strategy and demonstrate its efficacy on benchmark datasets without invoking ground truth.\n2. Given that only MarsGT appears to have undergone systematic tuning, while most baseline models run with default settings, could the authors provide a predefined and fully fair benchmarking protocol? Furthermore, tuning should be conducted under the same search space/stopping criteria using unlabeled rules. Label-based metrics (e.g., F1/NMI) should not be used for tuning.\n3. Regarding the EEL FISH dataset: The current explanation attributing discrepancies to \"over-clustering\" is insufficient. A concrete, step-by-step reconciliation analysis is required, including parameter sweeps and an identification of which preprocessing or resolution choices drive the divergence from original annotations.\n4. The framework currently appears to be an operationalization of classical graph statistics (clustering coefficients) rather than a fundamentally new algorithmic contribution. The authors should explicitly formalize the mathematical or conceptual novelty that distinguishes RareQ from deterministic propagation. Please provide a direct comparison against strong, generic baselines (e.g., standard label propagation variants, or Leiden/Louvain with resolution control) using the exact same kNN graph. This is necessary to prove that RareQ’s performance cannot be replicated by simple tuning of existing, widely-used pipelines.\n5. Rare population discovery is highly sensitive to technical noise. How robust is RareQ to upstream choices (e.g., number of PCs, HVG selection, batch correction methods, or k )? A systematic analysis is needed to show if/when RareQ produces spurious clusters due to density gradients or batch effects. The authors should include \"negative control\" experiments (e.g., permuted data, batch-only simulations, or datasets without true rare types) to quantify the False Positive Rate (FPR) of the method. Practical diagnostics for users to identify potential failure modes in unlabeled settings should also be provided."}, "author_response": "We thank the reviewer for the careful assessment of our revised manuscript. In this revision, we have (i) revised the hyperparameter discussion by removing potentially misleading statements and adding a label-free parameter selection/diagnostic framework, and (ii) expanded the analyses and comparisons to more explicitly demonstrate RareQ’s conceptual and empirical distinction from standard graph-clustering approaches. Below, we address each concern in detail.\n1. If RareQ is intended for unsupervised discovery, the reliance on label-dependent metrics (e.g., F1-score, NMI, Precision/Recall) to justify hyperparameter choices and \"high-accuracy regions\" is problematic. These metrics are unavailable in real-world discovery scenarios. Specifically, the suggestion in Supplementary Note 2 that users \"modify parameters accordingly\" to improve accuracy is impractical without ground-truth labels. The authors must provide a fully label-free parameter selection strategy and demonstrate its efficacy on benchmark datasets without invoking ground truth.\nResponse: Thank you for this insightful critique. Our original discussion in Supplementary Note 2 was intended to illustrate RareQ’s robustness across a broad parameter range; however, we acknowledge that the phrase “modify parameters accordingly” could be interpreted as requiring access to ground-truth labels, which are typically unavailable in real-world discovery settings. We have therefore removed this statement and implemented a fully label-free parameter selection framework that does not rely on external annotations. This framework is guided by three internal diagnostic criteria:\n(1) Rare-cluster credibility (Q): assesses local connectivity and internal cohesiveness.\n(2) Cluster independence (Qc): measures cluster-level separability by comparing intra-cluster connectivity to total connectivity.\n(3) Cross-parameter consensus (NMI ): represents the median pairwise NMI stability between the clustering obtained at a given parameter setting and those obtained under all other settings in the tested parameter grid, providing a label-free consensus measure across the parameter landscape.\nTo help users navigate the trade-offs among these diagnostics, we provide a ternary plot (Fig. R1a) summarizing their joint behavior across parameter combinations. Overall, RareQ shows robust performance across a broad parameter range. We further highlight the high-performing region and recommend selecting parameters from the Pareto-optimal frontier—settings that jointly maximize cohesiveness, independence, and stability.\nImportantly, parameters are chosen solely based on these label-free diagnostics. We then report F /NMI only as a post hoc benchmarking to validate the effectiveness of our proposed strategy. As shown in Fig. R1b, parameter settings within the trade-off region yield better classification of rare cells and major cell populations compared with settings from other regions of the ternary plot.\nFig. R1. Label-free parameter selection strategy for RareQ based on ternary plots. (a) The first column presents ternary plots summarizing the relative contributions of three label-free diagnostics (Q, Q , and NMI ) across different parameter combinations. c stability The second column shows zoomed-in ternary plots based on min–max scaling metrics. (b) Boxplots comparing rare-cell detection performance (F score) and global clustering accuracy (NMI) between parameter settings located near the center of the ternary plot and those farther from the center.\n2. Given that only MarsGT appears to have undergone systematic tuning, while most baseline models run with default settings, could the authors provide a predefined and fully fair benchmarking protocol? Furthermore, tuning should be conducted under the same search space/stopping criteria using unlabeled rules. Label-based metrics (e.g., F1/NMI) should not be used for tuning.\nResponse: We thank the reviewer for this comment. We first clarify that our benchmarking does not tune any method using ground-truth labels. All baselines, including MarsGT, were run with author-recommended default settings as described in the original publications and/or official documentation. This is a widely accepted benchmarking practice because it evaluates each method under configurations intended and validated by its developers rather than under arbitrarily chosen parameterizations. Ground-truth annotations are used only for post hoc evaluation, which is standard for benchmarking clustering and rare-cell discovery methods. To ensure full reproducibility and transparency, we had provided a predefined benchmarking protocol (introduced in the previous revision): identical preprocessing and input graphs where applicable, together with scripts and configuration files.\nRegarding the request to tune all methods under the same search space and stopping criteria, we respectfully note that a uniform search budget is not well-defined and not necessarily fair across heterogeneous algorithms. Their hyperparameters govern fundamentally different mechanisms (e.g., graph resolution for Leiden/Louvain, sketching and outlier scoring for FiRE, Gini cutoffs for GiniClust2, outlier thresholds for RaceID, and training dynamics for learning-based models), and their computational costs differ substantially. Imposing an identical search space/stopping rule can therefore bias the comparison and, in many cases, would replicate extensive sensitivity analyses already reported in the original method papers.\nThat said, to directly address the reviewer’s concern, we performed label-free parameters tuning using a predefined search grid and identical stopping criteria across methods (as similar to Fig. R1a). Parameter selection was guided only by our label-free diagnostics (i.e., internal cohesiveness/credibility, cluster independence, and cross-parameter consensus), with F /NMI reported strictly as retrospective evaluation. As shown in Fig. R2a, even under this standardized, label-free tuning protocol—allowing each baseline to select its best trade-off region—RareQ consistently achieves the strongest performance, consistent with our main conclusions.\nIn addition, we conducted a complementary sanity-check analysis to ensure our conclusions do not depend on any particular “best-case” configuration. Specifically, within each method’s author-recommended parameter range, we evaluated multiple representative parameter combinations (i.e., sampling diverse settings to span the plausible configuration space) and compared performance across methods. As shown in Fig. R2b&c, RareQ again consistently achieves the best results across these diverse parameter settings.\nTogether, both the requested label-free selection experiment and the broader parameter-range sanity check independently support the same conclusion: RareQ’s performance advantage is robust and not driven by asymmetric tuning or a single favorable parameter choice.\nFig. R2. Benchmarking RareQ against existing rare-population detection methods using a label-free parameter selection strategy across 20 scRNA-seq datasets. (a) Comparison of rare-cell detection performance (F score, precision, and recall) and global clustering accuracy (NMI) for the indicated methods, with parameter combinations selected using the ternary plot–based, label-free strategy. (b) Boxplots summarizing rare-cell detection performance across multiple sampled parameter combinations for each method. Each dot represents the result for a given method and parameter setting, with metric values aggregated across all rare cell types in the 20 datasets. (c) Same as (b), but showing global clustering accuracy measured by NMI.\n3. Regarding the EEL FISH dataset: The current explanation attributing discrepancies to \"over-clustering\" is insufficient. A concrete, step-by-step reconciliation analysis is required, including parameter sweeps and an identification of which preprocessing or resolution choices drive the divergence from original annotations.\nResponse: We thank the reviewer for this comment. Upon re-examining the original EEL FISH study, we realized that our previous interpretation of the reported '187 clusters' was inaccurate. We now understand that these clusters were derived from a cross-modal integration process (using Bonefight to map external scRNA-seq references, PMID: 36138169), which is conceptually distinct from a purely de novo clustering approach. Consequently, we withdraw the statement regarding over-clustering in the EEL study.\nWith this corrected understanding, it is clear that the previously perceived 'granularity discrepancy' does not stem from preprocessing or resolution choices, but from a fundamental difference in data analytical objectives. While the original study leveraged external reference data to achieve high-resolution mapping, our analysis demonstrates that RareQ can perform autonomous de novo discovery using only the measured 440-gene EEL FISH slice.\nTherefore, we focus on the fact that—even without an external reference—RareQ robustly recovers key hippocampal anatomical structures. This highlights RareQ’s ability to extract biologically meaningful spatial organization and identify rare-cell anchors from a targeted gene panel alone.\n4. The framework currently appears to be an operationalization of classical graph statistics (clustering coefficients) rather than a fundamentally new algorithmic contribution. The authors should explicitly formalize the mathematical or conceptual novelty that distinguishes RareQ from deterministic propagation. Please provide a direct comparison against strong, generic baselines (e.g., standard label propagation variants, or Leiden/Louvain with resolution control) using the exact same kNN graph. This is necessary to prove that RareQ’s performance cannot be replicated by simple tuning of existing, widely-used pipelines.\nResponse: As explained in our manuscript, the novelty of RareQ lies in characterizing a distinctive topological signature of rare populations on the single-cell kNN graph. Specifically, rare cells often form locally cohesive, clique-like microstructures with reduced external connectivity compared with major cell populations (Fig. 1b). RareQ operationalizes this topology into an explicit inference mechanism for rare-population discovery. To quantify this signature, we introduced a cell-wise neighborhood connectivity metric, Q.\nConceptually, Q is motivated by the same intuition as the clustering coefficient in small-world networks, i.e., measuring the “cliquishness” of a local neighborhood. However, Q is not a direct reuse of the classical clustering coefficient; rather, it is formulated to reflect the topology typically observed in single-cell kNN graphs, where rare populations are frequently embedded within dominant groups. Specifically, the classical clustering coefficient mainly measures how densely a node’s neighbors connect to each other (local compactness). In contrast, Q captures not only within-neighborhood cohesion but also the relative separation of that neighborhood from its surroundings: a cell attains a higher Q value when its neighbors are densely interconnected while exhibiting comparatively fewer connections to cells outside the neighborhood. This enables Q to jointly characterize local cohesion and relative independence (boundary sharpness), which is particularly informative for deeply embedded rare populations. Because the standard clustering coefficient is defined on undirected graphs, directly substituting it for Q on directed kNN graphs would be mathematically inconsistent. Q can be viewed as a kNN-graph–compatible reformulation of neighborhood “cliquishness” used by RareQ as an operational signal in inference rather than a post hoc statistic.\nImportantly, in RareQ, Q is not used merely as a passive descriptive statistic. Instead, it is embedded into the inference procedure. During Q-guided label propagation, high-Q local cliques act as stable topological anchors, directing label diffusion along highly consistent neighborhood structures until convergence. This propagation dynamics is mechanistically distinct from standard majority-vote label propagation or deterministic propagation variants. We further introduce a Q -guided recursive merging strategy to recover appropriate granularity for major populations without sacrificing rare-population separability, thereby preventing global over-clustering.\nTo directly address the reviewer’s question of whether RareQ can be replicated by simple tuning of widely used pipelines, we performed a same-graph benchmark against a strong generic baseline: Seurat’s default Louvain clustering. Using exactly the same kNN graph, we ran RareQ with default parameters and swept Louvain’s resolution across a broad range (0.5, 0.8, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 15, 20), evaluating rare-population detection (F , precision, recall) and global clustering consistency (NMI). We note that resolutions as high as 15–20 are far beyond typical Seurat usage and are included here solely as a stress test to examine whether aggressive tuning can reproduce RareQ’s behavior.\nThe results show that Louvain approaches RareQ’s rare-cell detection performance only at extremely high resolutions (15–20) (Fig. R3a). However, these settings induce pronounced over-clustering, excessively fragmenting major populations and substantially reducing global clustering consistency (markedly lower NMI) (Fig. R3a, b). In contrast, RareQ achieves both strong rare-population recovery and high global coherence without requiring extreme parameter settings. These results indicate that, while modularity-based clustering may theoretically approximate rare-population recovery through aggressive resolution tuning, it does so at the cost of generating many over-partitioned subclusters—making it difficult in practice to reliably distinguish true rare populations from artifacts of over-splitting in an unsupervised setting.\nIn summary, we emphasize that RareQ’s methodological novelty does not lie in proposing an entirely new notion of “neighborhood compactness,” but in:\n(1) Formalizing a cell-wise topological measure (Q) tailored to single-cell kNN graphs to capture rare-population cohesion and boundary sharpness, rather than directly adopting classical clustering coefficients;\n(2) Operationalizing this topology into a Q-guided, topology-aware label propagation mechanism in which high-Q cliques serve as stable anchors for diffusion and convergence; and\n(3) Coupling propagation with a Q -guided recursive merging strategy to suppress global over-clustering while preserving rare-population separability.\nConsequently, the characteristic behavior of RareQ cannot be reproduced by simple parameter tuning of existing modularity-based clustering methods.\nFig. R3. Performance comparison between RareQ and louvain clustering. (a) Comparison of RareQ and Seurat’s default Louvain clustering on identical kNN graphs. RareQ was applied using default parameters, while Louvain clustering was evaluated across a range of resolution values (0.5–20). (b) The number of detected rare cell clusters by RareQ and Louvain clustering.\n5. Rare population discovery is highly sensitive to technical noise. How robust is RareQ to upstream choices (e.g., number of PCs, HVG selection, batch correction methods, or k)? A systematic analysis is needed to show if/when RareQ produces spurious clusters due to density gradients or batch effects. The authors should include \"negative control\" experiments (e.g., permuted data, batch-only simulations, or datasets without true rare types) to quantify the False Positive Rate (FPR) of the method. Practical diagnostics for users to identify potential failure modes in unlabeled settings should also be provided.\nResponse: We thank the reviewer for this comment. In response, we conducted a systematic robustness evaluation across real datasets and designed dedicated negative-control simulations to quantify false positives. We also provide practical, label-free diagnostics to help users identify potential failure modes in unlabeled settings.\n(1) Robustness to upstream choices on real scRNA-seq datasets\nTo assess sensitivity to upstream analytical choices, we evaluated RareQ on 20 real scRNA-seq datasets while varying key preprocessing settings, including the number of PCs, HVG selection, and batch-correction methods. Overall, we observed that RareQ remains stable and achieves consistently strong performance when using sufficiently informative feature spaces, specifically when PCs ≥ 40 and HVGs ≥ 1,500 (Fig. R4a, b). Across these settings, RareQ’s rare-population detection and global clustering performance showed minimal variation, indicating robust behavior under commonly used preprocessing regimes.\nFor the k parameter (denoted as k.param to distinguish from RareQ’s neighborhood size k parameter) in kNN graph construction, we previously performed a dedicated sensitivity analysis in the first round of revision (Supplementary Fig. 24) by evaluating k.param in the range of 10–30. RareQ exhibited remarkable robustness across this range, with consistently stable rare-cell detection and overall clustering performance. Together with the simulation results presented in the manuscript and the additional analyses reported in Response 1, these evaluations support that RareQ is robust to typical parameter variations and maintains strong tolerance to technical noise.\n(2) Batch-effect stress test\nTo explicitly test whether RareQ may spuriously identify “rare populations” driven purely by batch composition, we designed a controlled simulation using the Splatter package (PMID: 28899397). We simulated datasets split into two batches and two biological cell types, with ~1,500 cells in total. We constructed a challenging scenario in which one batch within cell type 1 constitutes only 1% of all cells, whereas the remaining batch of the same cell type and both batches of cell type 2 each comprise ~33% of cells (Fig. R4c). We further varied the batch-effect strength by increasing Splatter’s batch.facLoc parameter from 0.1 to 1.0.\nWe then tested RareQ under (i) no batch correction and (ii) several commonly used batch-correction strategies, including CCA, RPCA, SCT, and Harmony (Fig. R4d). We asked whether RareQ would incorrectly call the rare batch (the 1% subset) as a rare “cell population”. Without batch correction, RareQ indeed identified the rare batch subset as a rare population, yielding FPR = 1 in this constructed batch-only rarity setting. With batch correction, RareQ largely avoided this failure mode when the batch effect was mild to moderate, and did not spuriously call the rare batch subset as a rare population (Fig. R4e).\nNotably, when the batch effect became extremely strong, batch correction may over-collapse biological distinctions such that batches across cell types become confounded, which can lead to mis-grouping. Among the tested approaches, CCA and SCT integration showed the most favorable behavior in preserving biological structure while mitigating batch-only artifacts. Overall, these results support a clear practical conclusion: when batch effects are present, applying batch correction is critical to prevent spurious “rare populations” that are driven primarily by batch composition rather than biology.\n(3) Negative-control dataset without true rare populations and FPR estimation\nTo further quantify the false positive rate under a “no true rare types” scenario, we constructed a negative-control dataset by extracting low-heterogeneity major populations from the UUO kidney dataset. This control dataset is designed to lack genuine rare populations. Across a range of PCs and HVG settings, RareQ produced near-zero FPR, indicating that RareQ does not tend to invent rare clusters in relatively homogeneous data when true rare populations are absent (Fig. R4f).\n(4) Practical diagnostics for unlabeled settings\nFinally, we provide practical diagnostics that users can apply when ground-truth labels are unavailable. First, because RareQ is topology-driven, the distribution of Q values offers an immediate label-free sanity check. In the UUO-kidney–derived negative-control dataset, Q-value distributions are similar across clusters, consistent with the absence of a topologically distinct rare clique (Fig. R4g). In contrast, when we spiked in a known rare population (e.g., JGA cells), the rare population exhibits substantially elevated Q values compared to other clusters, consistent with the intended topological signature that RareQ leverages (Fig. R4h).\nSecond, users can further assess biological plausibility using standard downstream analyses, such as differential expression/marker discovery (e.g., FindMarkers), to evaluate whether the detected rare cluster shows coherent marker programs consistent with a plausible biological state rather than technical artifacts.\nWe have incorporated these robustness analyses, negative controls, and user-facing diagnostics into the revised Supplementary information (see Supplementary Notes 4, 6, and 7), addressing the reviewer’s concerns about technical confounding and failure-mode detection in unlabeled settings.\nFig. R4. Robustness to upstream choices and batch-effect correction. (a-b) Robustness of RareQ to principle components (PCs) (a) and highly variable genes (HVGs) choices across real scRNA-seq datasets (b). (c) A simulation data with a batch-only rare cell cluster generated by the Splatter package. (d) False positive rate (FPR) assessment of RareQ on simulated data before batch correction and after batch correction with different methods across increasing batch-effect strengths; batch-confounded clusters are highlighted. (e) Cell populations identified by RareQ without batch correction and after batch correction using the CCA-based integration. (f) A negative-control dataset constructed from low-heterogeneity major cell populations extracted from the UUO kidney dataset; across a range of PCs and HVG settings, RareQ maintains near-zero FPR. (g) In the UUO-kidney–derived negative-control dataset, Q-value distributions are comparable across clusters. (h) Q-value–based diagnostic after spiking in a true rare population (JGA cells)."}]

Decision letter