(C) PLOS One
This story was originally published by PLOS One and is unaltered.
 . . . . . . . . . .



COEXIST: Coordinated single-cell integration of serial multiplexed tissue images [1]

['Robert T. Heussner', 'Department Of Biomedical Engineering', 'Computational Biology Program', 'Oregon Health', 'Science University', 'Ohsu', 'Portland', 'Oregon', 'United States Of America', 'Cameron F. Watson']

Date: 2025-09 

Multiplexed tissue imaging (MTI) and other spatial profiling technologies commonly utilize serial tissue sectioning to comprehensively profile samples by imaging each section with unique biomarker panels or assays. The dependence on serial sections is attributed to technological limitations of MTI panel size or incompatible multi-assay protocols. Although image registration can align serially sectioned MTIs, integration at the single-cell level poses a challenge due to inherent biological heterogeneity. Existing computational methods overlook both cell population heterogeneity across modalities and spatial information, which are critical for effectively completing this task. To address this problem, we first use Monte-Carlo simulations to estimate the overlap between serial 5μm-thick sections. We then introduce COEXIST, a novel algorithm that synergistically combines shared molecular profiles with spatial information to seamlessly integrate serial sections at the single-cell level. We demonstrate COEXIST necessity and performance across several applications. These include combining MTI panels for improved spatial single-cell profiling, rectification of miscalled cell phenotypes using a single MTI panel, and the comparison of MTI platforms at single-cell resolution. COEXIST not only elevates MTI platform validation but also overcomes the constraints of MTI’s panel size and the limitation of full nuclei on a single slide, capturing more intact nuclei in consecutive sections and thus enabling deeper profiling of cell lineages and functional states.

Multiplex tissue imaging (MTI) allows researchers to study tissue samples by measuring many proteins in individual cells while preserving their spatial organization, providing critical insights into the tumor microenvironment — the complex system of cells and structures surrounding tumors that influences cancer progression and treatment. However, current MTI platforms are limited in the number of biomarkers that can be stained simultaneously on a single tissue section due to panel size restrictions or assay incompatibilities. To expand the number of proteins measured, researchers often apply different panels to consecutive thin tissue sections. While image registration can align these sections, integrating them at the single-cell level remains challenging due to biological heterogeneity and differences in cell composition across slices. Existing methods often fail to fully incorporate both molecular and spatial information. To address this, we developed COEXIST, a computational algorithm that combines shared molecular profiles with spatial information to match cells across serial sections, enabling more complete and accurate tissue profiling. COEXIST improves MTI utility, enhances single-cell resolution, and offers opportunities for integrating MTI with spatial transcriptomics for deeper biological insight.

Competing interests: I have read the journal's policy and the authors of this manuscript have the following competing interests. The authors declare the following competing interests: G.B.M. is a SAB member or Consultant: for Amphista, Astex, AstraZeneca, BlueDot, Chrysallis Biotechnology, Ellipses Pharma, GSK, ImmunoMET, Infinity, Ionis, Leapfrog Bio, Lilly, Medacorp, Nanostring, Nuvectis, PDX Pharmaceuticals, Qureator, Roche, Signalchem Lifesciences, Tarveda, Turbine, Zentalis Pharmaceuticals. G.B.M. has Stock/Options/Financial relationships with: Bluedot, Catena Pharmaceuticals, ImmunoMet, Nuvectis, SignalChem, Tarveda, and Turbine. G.B.M. has Licensed Technology: HRD assay to Myriad Genetics, DSP patents with Nanostring. G.B.M. has Sponsored research with AstraZeneca. The other authors declare no competing interests.

Funding: This work was carried out with major support from National Cancer Institute (NCI) Human Tumor Atlas Network (HTAN) Research Centers at OHSU (U2CCA233280 to G.B.M and Y.H.C). Y.H.C is supported by R01 CA253860, U01 294548 and Kuni Foundation Imagination Grants. The resources of the Exacloud high-performance computing environment developed jointly by OHSU and Intel and the technical support of the OHSU Advanced Computing Center are gratefully acknowledged. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Data Availability: Data Availability As part of this paper all images at full resolution, all derived image data (e.g. segmentation masks), and all cell count tables are publicly released via the NCI-recognized repository for Human Tumor Atlas Network (HTAN; https://humantumoratlas.org/ ) at Sage Synapse (associated Identifiers: HTAN TNP – TMA, OHSU_TMA1_004-XX, https://www.synapse.org/Synapse:syn17022193/wiki/596115 ) where XX represents TMA core ID. The tonsil image dataset used for the label propagation experiment is available at https://www.synapse.org/MCMICRO_images and are available through the Human Tumor Atlas Network. The CRC WSI dataset is available through the Human Tumor Atlas Network (associated identifiers: HTAN TNP – SARDANA, HTA13_1_6, https://www.synapse.org/Synapse:syn18434611/wiki/590458 ). To request access, please log in or create a free Synapse account at https://www.synapse.org , and follow the procedures for data access as described on the respective dataset page. Code Availability All software used in this manuscript is detailed in the article’s Methods section and its Supplemental Information. The associated scripts are freely available via GitHub as described at https://github.com/heussner/coexist .

To address this, we first design a Monte Carlo simulation to estimate the theoretical overlap of single cells across consecutive tissue sections ( Fig 1B ). We then examine how LSA can be used with cell tracking algorithms, highlighting the importance of an approach based on tracking. Next, we introduce COEXIST, which leverages spatial information and shared biomarker expression to pair individual cells across consecutive tissue sections ( Fig 1C ). We apply our method to the Human Tumor Atlas Network (HTAN) Trans-Network Project tissue microarray (TNP-TMA) breast cancer (BC) dataset, where we phenotype cells using a combined tumor and immune panel of 42 unique molecular features, nearly doubling the number of features compared to a single panel. Our algorithm allows for the full reconstruction and characterization of 2.5D tissue, making it highly versatile for integrating spatial multi-omics data to analyze emerging spatial assays [ 39 ].

Alternatively, computational multi-omic integration methods aim to pair single-cell datasets and can be extended to spatial assays [ 28 ]. These methods assume similar cell states across datasets and use techniques such as canonical correlation analysis (CCA) [ 29 ] or variational autoencoders [ 30 ] to jointly embed single-cell representations and match them across modalities [ 31 ]. Linear sum assignment (LSA) [ 32 , 33 ] is frequently used as a framework to achieve this matching by constructing a cost matrix of potential cell pairs and identifying a set of pairs of minimal overall cost [ 34 , 35 ]. For example, MaxFuse [ 36 ], a recent approach backboned by LSA, is specifically designed for weakly linked modalities with few shared features, including proteomic, transcriptomic, and epigenomic information at single-cell resolution on the same tissue section. It generates nearest neighbor graphs, identifies anchor cells based on linked feature correlation, and then iteratively refines matches with graph smoothing and CCA. However, current multi-omic integration strategies are not compatible with spatial assays, as they fail to include spatial information and incorrectly assume homogeneity within the population, particularly when dealing with serial sections. Alternatively, cell tracking algorithms [ 37 ], designed to track cell movements across video frames, can be adapted to treat serial images as sequential frames but do not leverage the rich molecular profiles provided by MTI. Three-dimensional cell segmentation methods, such as the 3D extension of the Cellpose algorithm, similarly do not consider molecular profiles of cells across adjacent sections, are limited by the need for very thin tissue sections (2–3μm) and impose strict overlap thresholds to consider a cell tracked [ 38 ]. Therefore, it is essential to develop a single-cell integration method that combines molecular characteristics with spatial context to effectively integrate multi-modal or multi-panel datasets obtained from serial tissue sections.

(A) Visualization of shared (green) and unshared (yellow/blue) cells across serial tissue sections. (B) Monte-Carlo simulation estimates approximately 62% of cells are present in serial sections. (C) Overview of single-cell integration and downstream analysis: Cell-feature tables of serial tissue sections with shared and unshared features are integrated with an algorithm that accounts for the spatial and protein expression information, identifying the shared (green) and unshared (yellow, blue) cell populations. The shared cells undergo phenotyping using the combined panel, with resulting cell type labels being extended to the remaining unshared cells for each section. This propagation of cell type labels to the unshared cells is achieved via k-NN.

Early methods for tracking nuclei across serial electron microscopy images rely on manually drawn contours or thresholding to track objects across sections [ 18 ]. While modern image registration algorithms can align serial MTIs by identifying shared landmarks [ 19 – 22 ], achieving perfect alignment at the single-cell level is impossible due to inherent biological variation ( Fig 1A ). These types of approaches either require exceptionally thin sections, do not incorporate the full molecular profiles of the cells, or do not perform unified single-cell level analysis [ 8 , 20 , 22 , 23 ]. Furthermore, our previous work reported natural variations between serial sections in the image-to-image translation task, especially when using consecutive H&E sections to infer immunofluorescence intensity [ 19 , 21 ]. Recent platforms such as RareCyte Orion [ 16 ] or studies [ 24 ] are prioritizing multimodal assays such as CyCIF and hematoxylin and eosin (H&E) or high-plex protein combined with whole transcriptome sequencing in the same section to address this problem. Beyond MTI, many spatially resolved single-cell assays utilize serial tissue sections due to technological constraints [ 15 ], including Nanostring CosMX [ 25 ], spatial proteogenomics [ 26 ], and others [ 27 ].

Multiplexed tissue imaging (MTI) provides spatially resolved molecular profiling of single cells in thin tissue sections [ 1 ]. Validating MTI platforms and panels is essential for reproducibility [ 2 ], and is often achieved by comparing their performances on serial tissue sections. This comparison assumes uniformity across the sections and is conducted either at the pixel level by marker intensity or at the population level by cell phenotyping [ 3 – 6 ]. However, tissue sections are inherently heterogeneous, driven by factors including the sectioning process, the complexity of the tissue structure, or the presence of rare cell populations [ 7 – 10 ]. A recent study applying cyclic immunofluorescence imaging (CyCIF) to 40μm-thick melanoma sections found more than 90% of the nuclei were not fully intact within 5μm subsections [ 11 ]. This finding challenges the prevailing homogeneity assumption in bulk-level cell populations across serial sections. Another study of 4μm thick histological tissue sections showed their mean pixel-wise Pearson correlation was less than 0.6, underscoring their biological heterogeneity [ 12 ]. Moreover, fundamental challenges like spectral overlap and tissue degradation limit MTI panels to 40–100 proteins per sample [ 1 ]. By contrast, non-spatial technologies like CITE-seq and scRNA-seq can measure over 200 epitopes [ 13 ] and 1000 unique molecular identifiers [ 14 ] per cell, respectively. Although applying different MTI panels or spatial assays on serial sections allows for deeper molecular profiling, it sacrifices paired single-cell measurements. A computational approach to harmonize these datasets at the single-cell level would mitigate the challenge of unpaired single-cell measurements, improve rigorous validation of MTI technologies, and overcome MTI’s panel size limitations [ 15 – 17 ].

Validating across MTI platforms is critical for ensuring scientific reproducibility in spatial biology. We demonstrated COEXIST’s effectiveness in comparing MTI platforms using serial sections from normal jejunum core B4 in the TNP-TMA BC dataset processed with CyCIF and multiplex immunohistochemistry (mIHC) [ 17 ]. We compared the section bulk-level population compositions and found a strong Pearson correlation (r = 0.8) between the platforms ( Fig 5A ). This result is expected and has been demonstrated by other groups [ 4 ]. We performed the population comparison exclusively on the subset of shared cells identified by COEXIST, observing a similar correlation (r = 0.85) in compositions ( Fig 5B ). However, recognizing the limitation of bulk-level comparisons, we assessed the platforms at the single-cell level and found a moderate correlation, with Spearman correlations ranging from 0.26 (CD11B) to 0.81 (PanCK) ( Fig 5C ). The lower Spearman correlations observed at the single-cell level compared to the population level can be attributed to the limited presence of positive populations for certain markers (i.e., rare cell types). Additionally, we noticed the impact of the nonlinearity introduced by the IHC enzyme reaction, which affects signal conversion to expression levels.

Next, we propagated the rich phenotype labels to the unshared cells in each section using k-NN ( Methods ), showing increases in cell type diversity for each panel ( Fig 4E ). The immune panel inherits increased cell type resolution within the epithelial and stromal compartments, while the tumor panel shows a large expansion in immune cell type resolution. The resulting spatial distribution of cell types in the two sections appeared very similar after label propagation ( Fig 4E ). To quantify this similarity, we determined the cell type composition of all cells utilizing the k-NN propagation. The propagated cell composition remained consistent across panels when compared to their original counterparts ( S5 Fig ). We also documented recurrent cellular neighborhoods around cells shared by the two panels. Our recurrent neighborhood (RCN) analysis showed consecutive sections varied in the abundance of certain neighborhoods, such as RCN 3, 4, 5, and 8 ( Fig 4F ). This highlights the weakness in assuming spatial homogeneity across the tissue block. Other neighborhoods are similarly abundant across the two panels. One such neighborhood, RCN 7, is composed primarily of CD163 + Macrophages, CD31 + Endothelial cells, and CD20 + B-cells ( Fig 4G ). These lineage markers are spread across the two panels, and therefore the co-localization of these cell types would have been missed without the combination of the two panels. This RCN may suggest that vascular regions play an important role in immune cell infiltration. Further, this example highlights the utility of panel combination to uncover potentially biologically meaningful functional tissue units.

To further investigate cell identity switching, we examined a population of 138 cells that were classified as αSMA+ stromal cells in the tumor panel but reclassified as macrophages when the immune panel was incorporated ( Fig 4D ). For this population of cells, lineage marker intensities were ranked against the remaining cells in the core using each panel. In the tumor panel, αSMA and Vimentin emerged as the highest-ranked markers. Meanwhile, in the immune panel, markers such as CD14, MHC II, CD11b, αSMA, CD45, CD163, and CD68 were among the most prominently ranked, driving the identification of the macrophage population. On the combined panel (depicted in purple for markers originating from the tumor panel and green for those from the immune panel in Fig 4D ), CD14 and MHCII were notably more differentially abundant compared to αSMA. These differences likely stem from the absence of myeloid-specific lineage markers in the tumor panel. Without these lineage-specific markers, fluorescent spillover from nearby cells or high autofluorescence in other channels makes the misclassification of cell types more likely ( S4 Fig ). The integration of both panels effectively resolved this discrepancy. Apart from these markers, profiles across the cells were generally similar.

We then conducted a thorough phenotypic analysis focusing on triple-negative BC (TNBC) core E6 with combined tumor and immune CyCIF panels comprising 42 biomarkers ( S1 Table ). We identified 13 distinct cellular phenotypes, unveiling spatial heterogeneity at an unprecedented level ( Fig 4B and 4C ). However, clustering cells based solely on either the immune or tumor panels failed to capture their full phenotypic diversity ( Fig 4B ). For instance, the immune panel lacked markers such as CK14, CK17, PCNA, CD31, and Vimentin, which restricted its ability to delineate different epithelial sub-states, vasculature, or cells that may be undergoing Epithelial-to-Mesenchymal transition. Meanwhile, the tumor panel could not differentiate between CD8+ or CD4 + T-cells and lacked markers entirely for B-cells and all myeloid lineage cells including macrophages and monocytes. The combined panel rectified these shortcomings and even revealed a distinct stroma cell population. Shared markers across the panels clustered closely ( Fig 4C ) with moderate Spearman correlations ( S3 Fig ).

(A) Alluvial plot of high-level cell phenotype assignments of the matched cells generated by COEXIST, using either the tumor panel, immune panel, or combined panel for 39 BC subtype TMA cores. (B) Alluvial plot of refined cell types for E6 phenotyped with either the immune panel, tumor panel, or combined panel. (C) Heatmap of Leiden clusters used to define refined cell types of the matched cells in triple negative BC core E6 based on the combined panel. (D) Alluvial plot (top) and biomarker rank plot (bottom) highlighting a subset of tracked cells (n = 138) that switched cell type identities between the panels. For the biomarker rank plot, lineage markers were ranked for these 138 cells against the rest of the cells in the core. Purple represents biomarkers derived from the tumor panel, while green represents immune panel biomarkers. (E) Visual comparison of the cell phenotypes mapped to their spatial locations in the TNBC core consecutive sections (columns) before (top row) and after (bottom row) integrating them with COEXIST. (F) Recurrent neighborhood analysis heatmap showing cell composition per neighborhood and each neighborhood’s frequency in both sections. (G) Case study of RCN 7 from (F) showing image channels corresponding to RCN 7’s markers in both slides (top row) and the presence of RCN 7 after applying COEXIST (bottom row).

We performed a comprehensive analysis of all BC TMA cores, using combined tumor and immune CyCIF panels with a total of 42 unique biomarkers ( S1 Table ). We observed several discrepancies in the uni-panel phenotypes compared to the combined panel ( Fig 4A ). Specifically, cells profiled with the immune panel showed a predominant immune cell population, whereas the tumor panel did not exhibit such dominance. However, the combined panel revealed a more balanced population between these two cases, highlighting COEXIST’s ability to refine phenotypes. BC subtype-specific differences contribute to certain batch effects, making an integrative analysis challenging without rigorous normalization [ 42 ]. Nevertheless, high-level phenotyping showed evident biases towards the dominant cell population being profiled for a given panel, (i.e., immune or tumor cells), which COEXIST effectively corrects.

The experiment was performed on two datasets: (1) a region of interest (ROI) from a healthy tonsil specimen (~70,000 cells) ( Fig 3B ), and (2) a whole-slide image (WSI) of colorectal cancer (CRC) tissue with adjacent normal tissue (~480,000 cells) ( Fig 3C ). Results showed that as the proportion of tracked cells increased towards the theoretical maximum of 100% tracking across consecutive sections, label propagation performance, measured by F1 score, improved for both the tonsil ROI and CRC WSI datasets ( Fig 3D ). For k-NN propagation across 1000 random panels with 60% cell tracking, the mean F1 scores were 0.78 (SD 0.04) and 0.76 (SD 0.03) for the tonsil ROI and CRC WSI datasets, respectively. These scores approached the maximum theoretical mean F1 scores of 0.82 (SD 0.03) and 0.81 (SD 0.02). Importantly, this trend was consistent across varying anchor cell proportions, despite high variations between panels within each anchor cell proportion. Furthermore, both highly performing and poorly performing panels improved as the anchor proportion increased ( Fig 3E ). These findings highlight that in addition to the proportion of cells tracked across consecutive sections, the specific biomarker composition of the sub-panels plays a critical role in label propagation performance. This underscores the importance of carefully selecting markers to maximize the accuracy and reliability of cell-type assignments across tissue sections.

(A) Overview of single-slide validation experiment. Ground-truth cell type clusters were generated on all cells and all biomarkers for a tissue section dataset. Across a range of proportions, the cells were split into ‘tracked’ (green) and ‘untracked’ (yellow) matrices. Within each tracking proportion, the matrices were subset to a random biomarker panel half the size of the original panel. A k-NN classifier was trained on the ‘tracked’ cells and used to predict which ground-truth clusters an ‘untracked’ cell belonged to. The predictions were randomly subsampled and compared against the ground-truth labels from the ‘untracked’ cells. (B) Ground-truth Leiden clusters (n = 13) for a healthy tonsil dataset. Cells are annotated by Leiden cluster in a 2D UMAP representation of the full biomarker panel (left) and a spatial scatterplot (right). (C) Ground-truth Leiden clusters (n = 17) for a colorectal cancer whole-slide image dataset. Cells are annotated by Leiden cluster in a 2D UMAP representation of the full biomarker panel (left) and a spatial scatterplot (right). (D) Boxplot showing F1-score for the predicted cell type labels at each anchor cell proportion across 1000 random marker panels in the experiment for both the tonsil ROI dataset (orange) and the colorectal cancer WSI dataset (blue). (E) Histograms for both the tonsil ROI dataset (orange) and the colorectal cancer WSI dataset (blue) showing the distribution of Pearson correlation coefficients for all 1000 random panels. The correlation coefficient was calculated to test the correlation between anchor cell proportion and F1 score for each panel separately.

To ensure robust validation of our approach, we use full panel data and synthetically generated two panels to validate COEXIST for refining cell types and label propagation ( Methods ). In this setup, the cell types identified in the full panel serve as the ground truth, enabling us to assess how many anchors are required and evaluate the effectiveness of label propagation ( Fig 3A ). This experiment was designed to simulate realistic experimental constraints, where only a fraction of cells can be consistently tracked across consecutive slides, and a reduced subset of markers (50% shared markers in the case of the TMA dataset) is available for k-Nearest Neighbors (k-NN) propagation. We tested cell-type label propagation by synthetically generating sub-panels with a random selection of half of the biomarker from the full dataset, simulating the necessity of propagating labels using overlapping markers to untracked cells on consecutive sections.

While LSA generates quality matches, empirically fine-tuning neighborhood size is still challenging and subjective. Cell tracking is better suited for this problem because it directly leverages spatial and morphological information and does not require a radius parameter. However, cell tracking does not account for molecular information and is susceptible to erroneously matching cells that overlap but are merely adjacent in the lateral dimension. Therefore, we introduce COEXIST ( Methods ), which combines cell tracking with the molecular feature distance metric used in LSA. By integrating spatial, morphological, and shared feature information, COEXIST provides a robust framework for improving matching accuracy. As we demonstrated the weakness of LSA, we benchmarked COEXIST to the baseline cell tracking algorithm and found it consistently improved the shared feature correlations of the matches, often by more than 5% ( Fig 2C ). Small decreases in shared feature correlation for specific markers in some cores occurred due to artifacts or tissue loss in the round the markers were stained and imaged.

We also explored the robustness of spatially constrained LSA by measuring its cycle consistency [ 41 ], which assesses the algorithm’s ability to maintain consistency in assigning matches when swapping the order of the reference and target MTIs. To do this, we matched the original N-tracked reference cells with all potential target cells. Then we did the opposite by matching the original N-tracked target cells with all potential reference cells. We measured cycle consistency by calculating the fraction of shared matches between both sets relative to N ( Fig 2B ). Our findings revealed a rapid decrease in cycle consistency beyond 5μm neighborhoods, dropping to below 50% at 10μm. This underscores that LSA, widely used as the framework for multimodal integration [ 36 ], can yield unstable results when the cost matrix is asymmetric, leading to order-dependent solutions that vary based on the choice of reference and target MTIs and the cell counts in each. If the cost matrix were symmetric, linear sum assignment would yield order-invariant matches. However, in practice, serial sections rarely contain identical numbers of cells, and excluding unmatched cells from one side introduces asymmetry that destabilizes the pairing. In contrast, tracking-based algorithms inherently preserve cycle consistency by design.

(A) Tracked anchors outperform random anchors in reasonably sized neighborhoods. Error bars indicate the standard deviation of five seeds with the constrained case representing bulk-level or LSA without spatial constraints, a variant of MaxFuse. (B) Cycle consistency of LSA rapidly diminishes with increasing neighborhood size, whereas cell tracking ensures cycle consistency. (C) Heatmap showing the match correlation percent improvement when enriching the tracking algorithm with protein expression information for each TMA core and shared feature, with eight shared markers serving as ground truth for matching and integration. These markers are present in consecutive tissue sections within the dataset.

Within neighborhoods smaller than 50μm, or approximately 6 nuclei, where the true pairs are expected to reside, we observed higher quality matches for the N-tracked anchors than the N-randomly selected anchors ( Fig 2A ). Therefore, cell tracking is an essential precursor for attaining matches that are both high in quality and spatially coherent. However, when extending this process to all M (>N) reference cells, the resulting matches exhibit lower quality due to M exceeding N. This trend remains consistent across other cores ( S2 Fig ). As we know only a fraction of M is physically shared across sections, applying LSA without anchors is not ideal. To illustrate the bulk-level matching without considering spatial information, we employ LSA without using anchors or neighborhood constraints, which we refer to as a variant of MaxFuse without smoothing. MaxFuse’s graph smoothing effectively reduces variation within clusters while enhancing differences between them, enabling cluster-level rather than true single-cell matching across modalities. However, this approach is not suitable when true single-cell matching is required. Therefore, we emphasize that our benchmarking is conducted against MaxFuse without smoothing. Interestingly, match quality slightly surpasses this bulk-level correlation at neighborhoods larger than 200μm in the previous unanchored neighborhood expansion scenario ( Fig 2A ). This occurs because LSA optimizes linked feature correlation at the individual cell pair level, whereas we assess match quality at the matched population level. While the sum of costs is lower in the bulk-level case, suggesting potentially higher-quality individual matches, it does not translate to higher overall match quality as we measure it.

LSA achieves higher match quality compared to cell tracking given the same reference cells. This is because LSA directly optimizes shared feature correlation, whereas cell tracking does not, and considers matches across arbitrary distances, while cell tracking cannot ( S1A Fig ). Therefore, to implement LSA we introduced spatial constraints to the Jonker-Volgenant (JV) [ 40 ] algorithm to restrict potential matchings to a local neighborhood of size r. This modification also mitigates local registration discrepancies, overcoming possible failures in image registration. To create a fair benchmark for LSA and cell tracking, we first identified N anchor cells from a reference MTI that comprised M (M > N) cells using the cell tracking algorithm. This ensured that the benchmark reference cells were associated with spatially coherent target candidates ( Methods ). We then used our spatially constrained LSA implementation to identify target matches for the N anchors, iteratively expanding the radius parameter. Additionally, we conducted LSA by randomly sampling N anchor cells from the reference slide to test the importance of cell tracking in choosing suitable anchor cells instead of using random anchor cells. We then repeated the radial expansion process to create sets of matches.

Many multi-omics integration methods operate with the assumption that cell populations are shared across serial sections, allowing for the matching of populations by linking cells within the same molecular states. However, our simulation revealed applying this assumption directly to spatial assays and solely focusing on optimizing linked molecular feature correlation via LSA, such as in MaxFuse [ 36 ], does not guarantee spatial coherence of the matches ( S1A Fig ). Alternatively, adopting cell tracking algorithms that rely on cell morphology and location for this problem can guarantee spatial coherence [ 37 ] ( S1B Fig ). In our subsequent analysis, we investigate the use of both approaches for the integration problem and propose a novel algorithm that overcomes their respective limitations. We frame the problem as finding a single-cell matching between a reference and target image, representing consecutive MTIs. We define match quality as the mean Spearman correlation of the shared molecular features of the paired cells.

Our Monte Carlo simulation ( Methods ) indicates that approximately 62% of cells are shared across consecutive sections in the TNP-TMA BC dataset, which employs 5μm-thick sections and has a mean nuclear diameter of 7.8μm (with a standard deviation of 2.4μm). To simplify the simulation, we assumed nuclei to be spheres with normally distributed diameters, drawing nuclear size statistics directly from each image. At each simulation iteration, the nucleus diameter was sampled from a normal distribution parameterized by the observed nuclear diameter distribution (mean 7.8μm, SD 2.4μm) across all tissue cores, thereby accounting for the prevalence of smaller nuclei in the tissue. While the purpose of the simulation was to provide a rough estimate of the fraction of nuclei shared across consecutive sections, it aligned with recent literature [ 11 ] and closely matched results from the cell tracking algorithm. Additionally, we repeated the simulation using the nuclear size statistics and section thickness from another study of 3D tissue CyCIF [ 11 ]. Our simulation yielded comparable results to the ground-truth measurements: 89.3% of nuclei were not fully intact in 5μm thick sections (compared to >90%) [ 11 ], while 71.3-77.6% remained fully intact within 30–40μm sections (compared to 60–80%) [ 11 ]. Thus, our simulation provides a reliable estimate for the fraction of cells shared across serial sections, highlighting that both nuclear size and section thickness independently affect the probability of a nucleus appearing in multiple sections. The degree of cell colocalization varies substantially with section thickness and cell size, as shown in Fig 1B , where each bar represents a different TMA core. These results motivate a new approach to allow rigorous validation and integration of serial sections at the single-cell level.

Discussion

While staining a single tissue with both panels would be ideal, it is often impractical due to tissue loss and signal degradation caused by multiple rounds of CyCIF. To address this, our method introduces an innovative cell tracking and propagation framework that integrates complementary antibody panels across sections. This approach expands the phenotypic space without requiring all markers to be present in a single section. Building on this, we recognize that prevalent multi-omic integration methods often do not pair the same cells across modalities; instead, they rely on matching cells with similar molecular states using limited shared information. However, our study highlights that characterizing molecular states based solely on shared marker information for matching is not sufficient to capture the true spatial relationships between cells, and further, spatial information is crucial for precise single-cell matching.

To combat this challenge, we developed COEXIST and demonstrated its utility for precise single-cell matching across serial sections. Using COEXIST, we resolved differences between two orthogonal MTI panels with limited marker overlap and establish a set of refined cellular phenotypes on breast cancer tissues from a Human Tumor Atlas Trans-Network Project (TNP) Tissue Microarray. Further, we successfully establish detailed cell phenotypes for one triple negative breast cancer core. While individual panels can cause misclassification of some cell populations across serial sections, due to the absence of key lineage markers, COEXIST rectified these discrepancies. Our neighborhood analysis revealed unique and potentially biologically significant spatial structures, such as a neighborhood enriched for vascular cells and a diverse complement of immune cells, that could be captured only by the combined panel. Additionally, COEXIST mitigated the problem of overcalling specific lineages that is common in panels designed for specific cell types like immune cells, fibroblasts, macrophages, and malignant epithelial cells, thereby contributing to more accurate cellular phenotyping across diverse tissues.

We show COEXIST has O(n) time complexity with respect to the number of input cells, making it highly scalable to large MTIs (S6 Fig). However, the preprocessing step of computing the pairwise correlation distances is O(n2), which can be improved by only computing local distances rather than all possible distances. It is important to note that employing COEXIST necessitates the use of two consecutive slides, which should be considered in the context of experimental design and cost management. Despite the cost challenge, applying COEXIST two consecutive slides allows for the joint analysis of about 60% of the nuclei that are intact across two 5μm sections. COEXIST’s key assumption is that the shared cell population is representative of the broader tissue population. However, this assumption may not hold in cases where the tissue sections are thicker than the average cell diameter, potentially leading to a disproportionate representation of certain cell populations. For example, the immune cells, which tend to be significantly smaller than malignant cells, could be especially affected depending on how the section is cut in a 2.5-dimension. In this study, we focused on two consecutive sections, but this approach can be expanded to more sections while keeping the reduced sharing of cells in mind. Nonetheless, we analyzed three consecutive sections simultaneously by treating the middle section as our reference and then using pairs of two consecutive sections with COEXIST. Interestingly, the matched cells generated by COEXIST show only moderate shared feature correlations, which is unexpected considering the direct linkage of protein marker expression through tissue processing and antibodies. This contrasts with multi-omics, where the linkage between features is weaker and studies exploring differences in paired multi-omics measurements often reveal discrepancies [43,44]. Several factors contribute to these observed discrepancies; First, many shared markers identify rare cell populations, leading to most cells being “negative” (Fig 5C). Second, compartment-level differences may exist within the same cell, as the cells are being measured in different regions across consecutive sections. Third, technical batch effects are expected between MTI platforms, for instance, mIHC uses enzymatic labels for molecule tagging, whereas CyCIF employs immunofluorescent labels [45]. Comparatively, CyCIF-CyCIF integration exhibits stronger linear relationships in the expression levels of shared markers (S7 Fig). These factors should be considered in future multi-omics integration studies.

In addition to requiring two consecutive sections, COEXIST assumes high-quality image registration and that spatial proximity and shared marker expression serve as reliable anchors for matching. To test COEXIST’s robustness to registration errors, we rotated and translated one section and re-applied COEXIST, measuring how many original cell pairs were recovered. We found rotations of <0.3 radians and translations of <4 pixels resulted in 80% or more recovery (S8 Fig). Users should ensure quality registration through visual inspection and cross-referencing the matching percent with the Monte Carlo simulation framework. We do not impose strict assumptions about the distribution of cell features or spatial organization, acknowledging that these characteristics can vary widely across tissue types. Despite this flexibility, we recognize that this assumption may introduce biases, particularly in complex tissues with highly heterogeneous cell distributions. We also acknowledge the limitation of using shared feature correlation as a metric for match quality. However, this approach aligns with other publications [12,35] for assessing match quality. Publicly available datasets suitable for validating COEXIST across diverse tissue types remain sparse. To address this, we experimented using two single-slide datasets: one from a healthy tonsil specimen and another from a colorectal cancer (CRC) tissue sample. These datasets were used to validate the propagation of cell type labels under varying degrees of cell tracking efficiency and across random biomarker panels. The two tissue types represent very different spatial organizations — tonsil tissue is highly spatially stratified by cell type, while CRC tissue is more spatially mixed. Notably, both tissue types showed improved label propagation performance as cell tracking efficiency increased. To establish an upper bound for label propagation performance, we simulated perfect cell tracking by using 100% of cells for k-NN training. In this scenario, label propagation errors could only be attributed to missing key biomarkers in reduced panels or shared markers across consecutive sections. This finding highlights the importance of marker selection, especially for future multimodal applications where feature correspondence across modalities may be limited.

While existing 3D segmentation tools such as Mesmer and Cellpose have been successfully applied to thin sections (2–3µm), their applicability is more limited in thicker sections (e.g., 5 µm) or suboptimally prepared tissue where nuclear preservation across sections may be incomplete. COEXIST addresses these challenges by balancing biomarker profiles with spatial proximity to assign matches, improving flexibility in heterogeneous tissues. These features make COEXIST particularly well-suited for retrospective or opportunistic 3D reconstruction where sample preparation was not optimized for volumetric analysis. Moreover, COEXIST can complement existing 3D segmentation approaches by refining tracking results or propagating cell type annotations based on molecular similarity, especially in regions where segmentation continuity is disrupted. Looking forward, 3D MTIs with high lateral resolution can be used for ground truth validation of COEXIST when they become available.

Alternative computational integration methods employ strategies such as dimension reduction or graph smoothing to reduce complexity in high-dimensional data. However, MTI panels are relatively compact compared to broader single-cell assays and are often designed to measure orthogonal features, such as the tumor and immune panels in the TNP-TMA BC dataset. Therefore, those strategies are not required for our applications in this paper. For tasks requiring high-dimensional multimodal integration, COEXIST is highly compatible with joint embedding approaches [30,36], making it adaptable to a wide range of spatially resolved datasets. COEXIST is designed to allow researchers to fully leverage MTI technologies like CyCIF and mIHC and improve MTI assay development. While MaxFuse demonstrated that cell matching is possible between modalities with weakly linked feature spaces, we have shown that these matches are prone to be spatially incoherent. Given its flexible framework and ability to be applied to any spatially resolved assay, future use cases for COEXIST could include spatially aware integration across modalities with weakly linked features, such as spatial proteomics and spatial transcriptomics.

[END]
---
[1] Url: https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1013325

Published and (C) by PLOS One
Content appears here under this condition or license: Creative Commons - Attribution BY 4.0.

via Magical.Fish Gopher News Feeds:
gopher://magical.fish/1/feeds/news/plosone/