CellLENS enables cross-domain information fusion for enhanced cell population delineation in single-cell spatial omics data

Abstract

Delineating cell populations is crucial for understanding immune function in health and disease. Spatial omics technologies offer insights by capturing three complementary domains: single-cell molecular biomarker expression, cellular spatial relationships and tissue architecture. However, current computational methods often fail to fully integrate these multidimensional data, particularly for immune cell populations and intrinsic functional states. We introduce Cell Local Environment and Neighborhood Scan (CellLENS), a self-supervised computational method that learns cellular representations by fusing information across three spatial omics domains (expression, neighborhood and image). CellLENS markedly enhances de novo discovery of biologically relevant immune cell populations at fine granularity by integrating individual cells' molecular profiles with their neighborhood context and tissue localization. By applying CellLENS to diverse spatial proteomic and transcriptomic datasets across multiple tissue types and disease settings, we uncover unique immune cell populations functionally stratified according to their spatial contexts. Our work demonstrates the power of multi-domain data integration in spatial omics to reveal insights into immune cell heterogeneity and tissue-specific functions.

Extended Data Fig. 1 ∣. Detailed illustration of LENS-CNN and LENS-GNN-duo model architectures in the CellLENS pipeline.

The simplified versions are presented in Fig. 1. Additional details on CellLENS model architectures are described in the Methods section. Related code is deposited in the GitHub repository with documentation (see the Code Availability section for details).

Extended Data Fig. 2 ∣. Refined B cell subpopulations discovered by CellLENS in a tonsil CODEX dataset.

(A) Metrics evaluation of clustering performances on CODEX human tonsil tissue. Embeddings/representations of cells, from 7 different methods, were used as input: CellSNP representation, feature (protein expression table), concact (protein expression + neighborhood composition table), CCA representation, MOFA+ representation, SpiceMix representation, and MUSE representation. Clusters (for the calculation of Silhouette score, CH, and DB index) were generated with Leiden clustering using the same parameters. (B) UMAP visualization of the embedding and Leiden clustering result (cc: CellLENS clusters; fc; feature-only clusters). Left panel: CellLENS embedding; Right panel: feature expression table. (C) Zoom-in of the UMAP visualization of B cells (green), germinal center B cells (purple), and replicating non-GC cells (red) on CellLENS embedding (left) and feature-only embedding (right). Cells colored according to cell populations identified using CellLENS clusters in both UMAPs. (D) Zoom-in of the UMAP visualization of the same cell types shown in C, but colored by cluster numbers from CellLENS embedding (left) and feature-only embedding (right). CellLENS successfully separated replicating non-GC cells from GC B cells (c10 and c8), where feature-only failed (c8 partial). (E) B cell and replication-related marker expression (z-normed) heatmap. Left panel: clusters from CellLENS; Right panel: clusters from feature-only. (F) Visualization of cell type locations: B cells (green), GC B cells (purple), and replicating non-GC cells (red). Based on the location, replicating non-GC cells should not be mixed with GC B cells, as shown in the clustering result from feature-only.

Extended Data Fig. 3 ∣. Refined rare immune subpopulations discovered by CellLENS in a lymph node Xenium dataset (5k gene panel).

(A) UMAP visualization of the embedding and Leiden clustering result (cc: CellLENS clusters; fc; feature-only clusters). Left panel: CellLENS embedding; Right panel: feature expression table. (B) Heatmaps of top 10 mRNA marker genes for each cluster. Clusters from feature-only representation with Leiden clustering. Genes were identified via function ‘FindALLMarker’ in R package ‘Seurat’. (C) Same heatmap setup but clusters from CellLENS representation with Leiden clustering.

Extended Data Fig. 4 ∣. Supplementary plots illustrating the application of CellLENS on the human HCC CosMx-SMI data in Fig. 5.

(A) UMAP visualization of the embedding and Leiden clustering result (cc: CellLENS clusters; fc; feature-only clusters). Left panel: CellLENS embedding; Right panel: feature expression table. (B) Heatmaps of top 10 mRNA marker genes for each cluster. Genes were identified via function ‘FindALLMarker’ in R package ‘Seurat’. Left panel: clusters from CellLENS representation; Right panel: clusters from feature-only representation. (C) Visualization of spatial locations in the HCC tissue of different Macrophage subpopulations identified by feature-only representation. In each plot, cells in a target cluster are colored in red, HCC tumor cells are colored in grey, and other cells and empty spaces are colored in black.

Extended Data Fig. 5 ∣. Evaluation of CellLENS performance stability on adjacent CODEX tonsil sections: correspondence of CellLENS clusters identified individually on two sections.

CODEX tonsil tissue from two adjacent sections were used to evaluate the robustness of CellLENS performance. Shown here are Leiden clustering results (resolution = 1) based on CellLENS embeddings. In each of the three columns, left are clusters identified from CellLENS embedding of slide 1, and right are corresponding clusters in slide 2 aligned to their counterparts in slide 1.

Extended Data Fig. 6 ∣. Comprehensive benchmarking with additional methods.

A total of 12 methods are benchmarked here: CellLENS (our method). SpaGCN: Designed for spot-level spatial modalities (to identify spatial domains). StLearn: Designed for spatial modalities (to identify cell populations). SEDR: Designed for spatial modalities (to identify cell populations). MUSE: Designed for single-cell spatial modalities (to identify cell populations). SpiceMix: Designed for spatial modalities (to identify cell populations). BANKSY: Designed for spatial modalities (to identify cell populations or spatial domains). CellCharter: Designed for spatial modalities (to identify spatial domains). MOFA+: Designed for general modalities. CCA: General statistical procedure with canonical correlation analysis. Concatenation: Direct concatenation between feature and location matrix. Feature-only profile: Conventional way of cell type identification. We applied all 12 methods to five datasets presented in our manuscript (CODEX spleen, Xenium tonsil, CODEX tonsil, CODEX cHL, and CosMx HCC) and evaluated them using four different metrics (See Methods for details). Here we aggregated all the results across metrics and datasets into one summary figure. Each subpanel represents a specific metric (for example, Modularity score). The Y-axis indicates the average ranking of a method across benchmarking conditions (for example, K clusters or resolution numbers). On the X-axis, methods are arranged by their average ranking across all four metrics, such that methods on the left perform the best overall. In the summary figure, the rankings were averaged across all five datasets.

Extended Data Fig. 7 ∣. Comparison between embeddings generated by CellLENS and SpaGCN on CODEX mouse spleen data.

While initially designed for spot-level spatial-omics data, SpaGCN could run on the whole CODEX mouse spleen data ( ~ 50k cells). We compared the cell type delineation ability between the embeddings from CellLENS and SpaGCN. The figure setup is the same as described in Fig. 2. SpaGCN failed to identify the various B cell subpopulations captured by CellLENS.

Extended Data Fig. 8 ∣. Comparison between embeddings generated by CellLENS and SpaGCN on CosMx HCC data.

While initially designed for spot-level spatial-omics data, SpaGCN could run on the whole CosMx-SMI human HCC data ( ~ 50k cells). We compared the cell type delineation ability between the embeddings from CellLENS and SpaGCN. The figure setup is the same as described in Fig. 5. SpaGCN failed to identify the various Macrophage subpopulations captured by CellLENS.

Extended Data Fig. 9 ∣. Loss quantification of the CellLENS model when using different image feature extraction processes.

(A) We compared the CellLENS model training losses, across three CellLENS variations: 1) Default CellLENS, where the imaging feature extraction part is done by training an Alex-Net like CNN encoder (supervised by local cell type neighborhood composition vector). 2) CellLENS with a pre-trained ResNet50, where the image features were directly extracted with the pre-trained ResNet 50 model, flattened, and reduced to a vector with 128 dimensions. This vector is swapped with the original image feature vector obtained from the retrained Alex-Net model, and the rest of the CellLENS training process remains the same; 3) CellLENS with a pre-trained ViT (transformer), using a similar process as the pre-trained ResNet50 in (2).(B) We compared the losses from three CellLENS variations: 1) CellLENS default with Alex-Net as described above. 2) CellLENS but swapping out the Alex-Net architecture with a ResNet50 architecture, and retraining its weights (initialized at pretrained weights). 3) CellLENS but swapping out the Alex-Net architecture with a ViT architecture, and retraining its weights (initialized at pretrained weight). The model loss was calculated the same as described in the Methods section paragraph ‘Information retrieval efficacy evaluation of the LENS-GNN duo module’. In these cases, we implemented a 80/20 train test data split. Retraining was only done on the train data, and loss values were calculated on test data.

Extended Data Fig. 10 ∣. Memory and run-time benchmarking for training CellLENS on datasets of different sizes.

We ran CellLENS on increasing numbers of cells: ~10k, ~50k, ~150k, ~500k, ~1.5 mil, and recorded the memory and run-time statistics. The benchmarking was performed on an NVIDIA A5000 GPU with an SSD disk (PCIe 4.0) for I/O. Left panel: GPU memory usage during CNN and GNN training stages of CellLENS, across variable cell numbers. Middle panel: Runtimes of CellLENS CNN training stage. Recorded run-time for both default LENS-CNN or LENS-CNN with a ViT architecture. Right panel: Run-time of CellLENS GNN training stage.

Fig. 1 ∣. Illustration of the CellLENS pipeline.

CellLENS is compatible with imaging-based spatial omics modalities with single-cell or finer resolutions (for example, CODEX, Xenium and CosMx). Information from three domains is extracted from each individual cell and its surroundings: (1) single-cell expression profile (for example, measured protein or mRNA features); (2) single-cell location information (for example, cellular neighborhood composition); and (3) singlecell local tissue image information (for example, local images from nuclear and membrane channels). CellLENS takes these three types of information as input. It first uses a CNN encoder to extract features from images of local tissues surrounding each cell. Next, two separate GNN models are constructed in parallel: (1) a ‘spatial GNN’, where each node represents a cell, with the initial node vector assigned as CNN-extracted local image features; nodes are connected according to spatial adjacency; and (2) an ‘expression GNN’, where each node represents a cell, with the initial node vector assigned as the expression profile; nodes are connected according to expression similarity. The two GNNs are connected by an overarching MLP head that combines the message-passing outputs of the two GNNs for predicting the target vector of each cell, that is, the concatenation of the cell’s feature-based population identity (one-hot) and its neighborhood composition (percentage) vectors. After training, the last layers of the two GNN models are extracted, combined and reduced (via singular value decomposition) to form the final, tri-domain integrated representation vector for each cell. This multi-domain fused representation vector is then used in the downstream analysis for cell-type identification purposes, which is compatible with commonly used unsupervised clustering methods (for example, Leiden clustering). Detailed illustration of the model architecture can be found in Extended Data Fig. 1. Cell illustrations created using BioRender.com.

Fig. 2 ∣. Refined B cell subpopulations discovered by CellLENS in a healthy mouse spleen CODEX dataset.

a, Metric-based evaluation of cell population delineation performances on CODEX mouse spleen tissue. Representations of cells, from seven different methods, were used as input: CellSNP; feature (protein expression table); concact (protein expression + neighborhood composition table); CCA; MOFA+; SPICEMIX; and MUSE (Methods). Five batches, each with 10,000 randomly selected cells, were tested. The solid line indicates the average and the shading indicates the 95% confidence interval (CI) of the scores. b, Uniform manifold approximation and projection (UMAP) visualization of representations and Leiden clustering (cc, CellLENS clusters; fc, feature-only clusters). The cell types of the CellLENS or feature-only clusters were annotated based on the average expression profiles of the clusters. Left, CellLENS representation. Right, feature expression. The dashed line indicates the B cell subpopulations. c, UMAP visualization and B cell-related protein expression profiles of clusters annotated as B cells. Left, B cell clusters from CellLENS and their expression heatmap. Right, B cell clusters from feature expression and their expression heatmap. d, Comparison of spatial locations of different B cell clusters identified using CellLENS representation versus feature expression in the spleen tissue. In each plot, the red dots indicate cells from a specific cluster; the green outlines indicate GC boundaries. Top, Spatial locations of B cell subpopulations identified by the CellLENS representation clusters. Bottom, Spatial locations of B cell subpopulations identified using the feature expression clusters. Annotations shown in parentheses summarize the key features and spatial characteristics of the respective clusters. NK, natural killer (cell).

Fig. 3 ∣. Refined immune cell subpopulations discovered by CellLENS in a tonsil and LN Xenium dataset.

a, Metrics evaluation of clustering performances on CODEX human tonsil tissue. Embeddings and representations of cells, from seven different methods, were used as input: CellLENS; feature (protein expression table); concact (protein expression + neighborhood composition table); CCA; MOFA+; SPICEMIX; and MUSE. Clusters (for the calculation of the silhouette score, and the Calinski–Harabasz and Davies–Bouldin indices) were generated with Leiden clustering using the same parameters. The solid line indicates the average and the shading indicates the 95% CI of the scores. b, Zoom-in of the UMAP visualization of B cells (green), GC B cells (purple) and replicating non-GC cells (red) on CellLENS embedding (left) and feature-only embedding (right). Cells are colored according to the cell populations identified using the CellLENS clusters in both UMAPs. c, Zoom-in of the UMAP visualization of the same cell types shown in c but colored according to the cluster numbers from CellLENS embedding (left) and feature-only embedding (right). d, Visualization of the B-proliferating cluster locations. Left, B-proliferating cells from the feature-only representation clustering (fc3), showing a mixed spatial pattern. Middle, B-proliferating cells from the CellLENS representation clustering (cc7), showing a GC spatial pattern. Right, B-proliferating cells from the CellLENS representation clustering (cc15), showing a non-GC spatial pattern. e, General cell-type locations on the Xenium human LN data (5k gene panel). f, Uniquely identified rare immune cell populations using CellLENS and their corresponding locations with the GC boundaries outlined in black.

Fig. 4 ∣. Refined T cell subpopulations in TMEs discovered using CellLENS in a cHL tumor CODEX dataset.

a, Metric-based evaluations of cell population delineation performances on CODEX human cHL tissue. Representations of cells, from seven different methods, were used as input: CellLENS; feature (protein expression table); concact (protein expression + neighborhood composition table); CCA; MOFA+; SPICEMIX; and MUSE (Methods). Five batches, each with 10,000 randomly selected cells were tested. The solid line indicates the average and the shading indicates the 95% CI of the scores. b, UMAP visualization of the representations and Leiden clustering (cc, CellLENS clusters; fc, feature-only clusters). The cell types of the CellLENS or feature-only clusters were annotated based on the average expression profiles of the clusters. Left, CellLENS representation. Right, feature expression. c, Visualization of cell-type spatial locations in the cHL tissue, colored according to the annotations on the CellLENS clusters. The black regions are empty spaces. The white outlines indicate the borders of the cHL tumor regions. d, Visualization of the spatial locations of different CD4 T cell subpopulations identified by the CellLENS representation clusters. The black outlines indicate the borders of the cHL tumor regions.

Fig. 5 ∣. CellLENS-enabled delineation of biologically distinct macrophage subpopulations in an HCC tumor CosMx SMI dataset.

a, Metric-based evaluation of cell population delineation performances on CosMx SMI human HCC tissue. Representations of cells, from seven different methods, were used as input: CellSNP; feature (protein expression table); concact (protein expression + neighborhood composition table); CCA; MOFA+; SPICEMIX; and MUSE (Methods). Five batches, each with 10,000 randomly selected cells were tested. The solid line indicates the average and the shading indicates the 95% CI of the scores. b, Visualization of spatial locations of different cell populations, including all cell types (first panel colored according to the cell-type annotation obtained from the CellLENS clusters; black regions indicate empty spaces) and different macrophage subpopulations identified by the CellLENS representation clusters (second to fourth panels). In each of the second to fourth panels, all tumor cells, and macrophage cells from a specific CellLENS cluster, are colored, while other cells and empty spaces are shown in black. c, Volcano plot of differentially expressed genes between the CellLENS cc6 cluster and other macrophage clusters. Statistics were generated using Limma. d, Comparison of module score values between CellLENS cc6 (n = 2,569) and all other macrophage cells (n = 5,089). ‘M1-like’ and ‘M2-like’ scores were calculated using genes from ref. . Splenic macrophage-specific ‘pro-inflammatory’ and ‘immunoregulatory’ scores were calculated using the genes from ref. . An unpaired two-sided Wilcoxon rank-sum test was used to determine the P values. The box plot indicates the 25th to 75th percentiles (box boundary), the median values (line in the box), and the minimum and maximum values (whiskers). e, Visualization of the spatial distribution of all macrophages and their respective LR interaction detection score levels. The detection score was calculated based on significant LR interaction pairs between macrophages and tumor cells. f, Top ten most frequent LR interaction pairs associated with CellLENS cc6 macrophages. g, GEP usage scores among tumor cells, stratified according to infiltration (by macrophage) status. Infiltrated tumors n = 5,953; other tumors n = 18,737. An unpaired two-sided Wilcoxon rank-sum test was used to determine the P values. The setup of the box plots is the same as in d. LSEC, liver sinusoidal endothelial cell.