EMR samples were collected during clinically indicated endoscopic procedures in patients referred to Cambridge University Hospitals NHS Trust (Addenbrooke’s Hospital, Cambridge, United Kingdom) for the treatment of BE-related neoplasia. The study was approved by the Institutional Ethics Committees, and all subjects gave individual informed consent for the use of their tissue samples for research purposes (REC 01/149). Briefly, after the endoscopic procedure, EMR specimens were pinned down on cardboard using 20 mm Agani needles and fixed overnight in formalin and processed for histology as per clinical standard. Approximately 5 μm thick FFPE sections were cut, and adjacent sections were stained for hematoxylin and eosin (H&E) for histopathological evaluation.

The H&E slides were scanned and printed out. Two experienced specialist gastrointestinal pathologists, MT and AM, independently assessed the H&E slides using the definitions and histological criteria for normal esophagus, Barrett’s metaplasia, the grades of Dys, and esophageal adenocarcinoma recommended by The Royal College of Pathologists (23) and guidelines by the British Society of Gastroenterology on the diagnosis and management of Barrett’s esophagus (24). There were no discrepancies between the two pathologists when reviewing the samples presented in this study. Indefinite Dys was excluded; high and low grades of Dyswere grouped together as Dys. All local pathological grades were marked on the corresponding H&E printout. ROIs were then determined based on the marked H&E printout. Unstained FFPE slides that were adjacent to the assessed H&E slides were selected to proceed for mIHC.

The staining was carried out on Bond RX platform (Leica System), following modified standard ‘IHC protocol F’ and Heat Induced Epitope Retrieval 1 (HIER1) for 20 min using the BOND Polymer Refine Detection kit. The ‘IHC protocol F’ was modified by 1) inserting three extra 5 min washing steps after HIER; 2) replacing the step ‘Mixed DAB Refine’ with ‘Bond Open Container’, which was supplied with freshly made 3-Amino-9-Ethylcarbazole (AEC); and 3) removing the step of ‘Hematoxylin’. One cycle of mIHC was carried out as one staining protocol, following the standard operation of the manufacturer’s protocol, except a) in the first cycle, the ‘Dewax’ step was selected, and b) in the last cycle, the ‘Hematoxylin’ step was added back as a sole staining step.

After each cycle of staining, slides were unloaded from the Bond machine and temporarily mounted in Tris-Buffered Saline with 0.1% Triton X-100 (TBST) and imaged using Zeiss AxioScan Z1 at 20x brightfield. After imaging, AEC was removed by ethanol by 2 brief washes in distilled water, 1 wash in 70% ethanol, and 1 wash in 100% ethanol for 3.5 min. The slides were then rehydrated through 2 min incubation in 70% ethanol, 1 min incubation in 30% ethanol, and then 4 washes in distilled water. The slides were then left to rest in TBST for the next cycle.

The raw images were generated in Carl Zeiss Image (CZI) format with resolution of 0.22 µm per pixel by Zeiss AxioScan Z1 at 20x brightfield. The raw CZI images were processed using the Zeiss software, Zen Lite, to generate ROI images in TIFF format according to the pathologist’s grading ( Figure 1B ). Each ROI has 11 raw Tag Image File Format (TIFF) images including nuclear staining and 10 biomarkers.

For each ROI, image coregistration was performed using a MATLAB-based script: ‘register_crop_batch4tif_CRUK_OHSU_mIHC.m’, which generated coregistered images shown in Figure S2B . The deconvolution of the 10 marker images (AEC stained) and the nuclear segmentation of the nuclear image were performed via the ImageJ macro: ‘AEC_CMYK_Seg_Overlay_Batch_CRUK_OHSU_mIHC.ijm’. This pipeline generated a simulated IHC image of each marker ( Figure S3A ) and monochrome images for simulated immunofluorescence ( Figure S3B ). It also generated the nuclear mask, a TIFF image with each single cell segmented, and a raw CellProfiler Output (CPOUT) file for FCS Express™ Image Cytometry (De Novo Software, Los Angeles, CA, USA) for manual gating, which were described in (15) and in Figure S4 . The cell numbers of all gates were exported as comma-separated values (CSV) files ( Figure 2B ) for the analysis of cell density and proportion and then used in the sparse subspace clustering (SSC) approach and spatial analyses.

In this approach, to identify the groups of cells with similar features, we facilitated the characterization of biologically significant cellular characteristics where representative data points are selected as landmarks, representing the original data points as linear combinations; X=[ x ₁, ⋯, x _N ] is an m × N data matrix where x_i represents a feature vector corresponding to the i-th cell, m represents the dimensionality of a feature vector x_i , and N is the number of segmented cells.

The main idea of SSC takes advantage of the self-expressiveness property of the data; for instance, each data point in a union of subspaces can be efficiently reconstructed by a combination of other points in the dataset. More precisely, each data point (or the expression of proteins for a segmented cell) x_i can be written as x_i = X c_i where c _i ≜[ c _i1 c _i2⋯c _N ] ^T and the constraint c_ii = 0 eliminates the trivial solution of writing a point as a linear combination of itself. In other words, the matrix of data points (i.e., the expression of proteins of all segmented cells) X is a self-expressive dictionary in which the expression of the proteins of an individual cell can be written as a linear combination of the expression of proteins of other cells. In general, as the representation of x_i in the dictionary X is not unique, among all solutions, we seek a sparse solution, c_i , whose nonzero entries correspond to the data points from the same subspace as x_i . By doing this, we identify a subspace-sparse representation where a sparse representation of a data point finds points from the same subspace where the number of the nonzero elements corresponds to the dimension of the underlying subspace. One can restrict the set of solutions by minimizing an objective function such as the l ₁-norm of the solution as follows:

After solving the proposed optimization program, we obtain a sparse representation for each data point whose nonzero elements ideally correspond to points from the same subspace. The next step of the algorithm is to infer the clustered group of the data into different subspaces using the sparse coefficients. The clustering of data into subspaces then follows by applying spectral clustering to the similarity graph in which the nodes that correspond to points from the same subspace are connected to each other and there are no edges between the nodes that correspond to points in different subspaces (k). SSC manages complexity by selecting a few representative points as landmarks so that the spectral embedding of the data can be efficiently computed with the landmark-based representation.

We used the Single-Cell Image Analysis Package (SCIMAP) open-source python library (https://github.com/labsyspharm/scimap) to quantify the average shortest distance between reference and target cells (scimap.tl.spatial_distance) and implemented neighborhood enrichment analysis to compute how likely cell types are found next to each other compared to a random background (scimap.tl.spatial_interaction). This uses a permutation test to compare the number of interactions between all cell types in a given image to that of a matched control containing randomized cell phenotypes. Thus, it enables the unbiased and systematic study of cell–cell interactions present in all the tissues of a sample cohort. By doing this, it determines the significance of cell–cell interactions, reveals enrichments (red) or depletions in cell–cell interactions (blue) that are indicative of cellular organization, and statistically non-significant results (gray). For parameters, we used the k-nearest neighbor algorithm (knn = 10) to identify the neighbors for every cell and the number of permutations with 1,000 as default. P-values are calculated by subtracting the permuted mean from the observed mean divided by the number of permutations as described in (25). For analyzing multiple images together, a cluster map shows the average cell–cell interaction across all images.

For Figures 4 , 5 , 7 , we used the Mann–Whitney–Wilcoxon test two-sided with Bonferroni correction. In a given cell population, we treated its density in the four different stages: normal esophagus with squamous epithelium (NSQ), non-dysplastic BE (NDBE), Dys, and EAC, as four independent groups. One sample point represented one ROI, and the p-value was calculated individually between two selected stages, e.g., NSQ vs. NDBE or Dys vs. EAC, whereby each stage contained a number of non-overlapping ROIs. For Figures S6–S8 , those cell subsets have high variations, whereby we did not particularly aim to investigate the change between each disease stage. We therefore chose the nonparametric Kruskal–Wallis one-way ANOVA test in R to determine if a given cell subset is statistically significantly different by at least one disease stage. For Figures 8 –F , we applied the Mann–Whitney-Wilcoxon test two-sided with Bonferroni correction for multiple comparison with the 5% significance level [open-source python library (26)].

The principles underlying mIHC have been previously described (22, 27). Briefly, each signal marker is detected by chromogenic-based (3-amino-9-ethylcarbazole, AEC) IHC in a sequential and multi-cycle manner, followed by digital imaging. After imaging, AEC is removed by ethanol and antibodies are stripped by heated citrate buffer before entering the next cycle ( Figure 1A ). Here, we modified the process to enable automated staining using the Leica BOND Automated IHC Stainer platform where each cycle could be performed as one standard staining program, and up to 30 slides could be stained in the same batch. We also confirmed complete removal of markers from previous cycles using the BOND standard heat-induced epitope retrieval step ( Figure S1 , Table 1 and Method).

We generated a series of digitized images for each marker representing all the EMR tissues, which were assessed independently by two pathologists for local disease grades. This allowed us to identify and acquire 40 ROIs spanning from an NSQ, NDBE, Dys, and EAC ( Figure 1B ; Table 2 ). Each ROI measured approximately from 0.5 × 1 mm to 1 × 2 mm, contained 11 single channel images that corresponded to the nuclear staining, and 10 biomarkers including CD45, CD3, CD8, FOXP3, CD20, CD68, CSF1R, CD163, CD1C, and KI67, respectively ( Figure 2A ). Previously, images were co-registered manually by identifying a “fiducial” point through all images. Here, we established an automated image registration algorithm based on the distinct tissue pixel gradients of the ROI and feature matching to coregister the images of different markers of the same ROI (21, 28). In the computational process, first, we extracted feature descriptors, matched features by using their descriptors, and retrieved the locations of corresponding points for each image. Then, we estimated transformation corresponding to the matching point pairs and recovered the scale and angle by using the geometric transformation. This algorithm enabled a fully automated coregistration and batch processing ( Figure S2 ). To assess the 10 biomarkers and their expression on every single cell, coregistered image stacks were processed using a watershed-based segmentation in FIJI (Fiji is Just ImageJ), which segmented single cells based on the hematoxylin nuclear staining, and then quantified the chromogenic intensities of AEC staining for each biomarker. Image analysis pipelines also generated multi-channel-merged IHC images that were reassessed by pathologists ( Figure S3A ), with monochrome images used to generate multi-channel images to visualize selected biomarkers ( Figure S3B ). Using the automated IHC staining and automated computational image analysis pipelines, we acquired raw data for 712,252 cells from the 40 ROIs with quantitative expression values from 10 biomarkers on a single-cell basis, including the cells’ shape–size measurement and spatial coordinates within the tissue ( Figure 2B ).

To identify key immune cell populations from the ROIs, we loaded raw data into FCS Express™ Image Cytometry (see Materials and Methods). Similar to flow cytometry, each data point corresponds to a cell projected in a biaxial plot that allows for quantitative assessment and selection (gating) based on the markers of choice. In addition, spatial gates could be directly annotated on the tissue images to select cells within a specific region. We therefore specifically focused on the BE/EAC infiltrating immune cells via gating cells within the BE or EAC epithelium and adjacent 200 µm wide stromal margin in all ROIs ( Figure S4 ). We then gated and assessed 17 key cell populations of lymphoid and myeloid lineages by adapting the previously reported gating strategies (20, 27) ( Figure S4 ; Table 3 ), including but not limited to CD45⁺ total immune cells, CD45⁺CD3⁺ T cells, CD45⁺CD3⁺CD8^-FOXP3⁺ T regulatory (Treg) cells, CD45⁺CD3^-CD20⁺ B cells, CD45⁺CD3^-CD20^-CD1C^-CD68⁺CSF1R⁺ myelomonocytic subsets, and the reported myeloid dendritic cell populations (13, 17, 18), marked by CD45⁺CD3^-CD20^-CD1C⁺ (hereinafter referred to as the CD1C⁺ subset). Gated cells were then visually overlaid onto original IHC images to confirm gating accuracy ( Figure S5 ).

In image cytometry gating profiles, we observed that the total cells were shifting toward higher CD45⁺ density with disease progression; there was also an obviously increased presence of CD45⁺CD3⁺CD8^-FOXP3⁺ Treg cells and the CD45⁺CD3^-CD20^-CD1C^-CD68⁺CSF1R⁺CD163⁺ myelomonocytic subset (hereinafter referred to as CD163⁺ myelomonocytic cells) ( Figure 3 ). We then plotted the cell densities (cell number per square millimeter) of all populations according to the disease grade of each ROI ( Figures 4 and Figure S6 ). Among all, the total CD45⁺ immune cells exhibited a clear increasing trend as the disease advanced ( Figure 4A ), which was in line with the overall pro-inflammatory microenvironment of BE/EAC development (29). Treg cell density showed a steady increase from NSQ to EAC ( Figure 4B ), where two significant increases were observed in NSQ (approx. 30 cells/mm²) versus NDBE (approx. 115 cells/mm²) and from Dys (approx. 76 cells/mm²) versus EAC (approx. 160 cells/mm²). A similar increasing trend was also observed in CD163⁺ myelomonocytic cell subsets, which reached the highest level in Dys (580 cells/mm²) and appeared to remain constant in EAC ( Figure 4C ). The overall trend in the CD CD163^– myelomonocytic cell subset was less clear, whereby it had significantly higher density in all disease grades compared with NSQ but appeared to reach the highest level in NDBE ( Figure 4D ). Interestingly, the abundance of both total T cells (CD45⁺CD3⁺, Figure 4E ) and CD8⁺ T cells (CD45⁺CD3⁺CD8⁺, Figure 4F ) did not show a clear trend, with similar observations made in the CD1C⁺ subset, too ( Figure 4G ).

In addition to cell density, the immune complexity is also reflected by cell composition; we therefore assessed for each cell group ( Figures 5 , S7 ). Notably, the proportion of CD8⁺ T cells decreased dramatically as the disease advanced ( Figure 5A ). In addition, the proportion of CD163⁺ myelomonocytic cell subset in total myelomonocytic cells significantly increased, from approx. 15% in normal NSQ and NDBE to approx. 50% in Dys and EAC ( Figure 5B ). The proportion of total CD45⁺ immune cells ( Figure 5C ) and Treg cells ( Figure S7 ) were in line with the cell densities that had a significantly increasing trend, but no clear trend was observed in the CD1C⁺ subset ( Figure 5D ).

The hierarchical gating strategy in image cytometry was useful in identifying immune cell populations from our multiplex image dataset. However, it relied on a priori information, such as canonical or known marker combinations and expertise in gating and confirming the cell classifications for each ROI ( Figures S4, S5 ). In addition, because the gating was performed manually on biaxial plots, the data were only assessed sequentially for up to two markers at a time; therefore, in high-dimensional data, such as this dataset, other important markers that describe cell phenotypes may be missed. In the case of using multiple biomarkers to evaluate more than 10 parameters, manual analysis via this gating strategy becomes a significant expenditure of time. For example, the interrogation of ₁₀C ₂ = 45 biaxial plots is required for the evaluation of 10 biomarkers. In addition, while robust, the inherent subjectivity of manual gating diminished its practical application for the clinical evaluation of specimens. Therefore, efficient and objective interpretation of mIHC-stained image entails several challenges limiting broad methodologic and clinical dissemination.

To overcome the limitation of the manual gating strategy and handle both biological heterogeneity (e.g., various cell types) and high redundancy in feature representation, we adopted the SSC approach by extending a similar concept of sparse coding used in our previous work (22) (see Materials and Methods). The proposed SSC approach enables objective and automated cell clustering via a simultaneous assessment of all 10 markers from all ROIs. In addition, the SSC approach has a better overall performance than other approaches such as principal component analysis (PCA) when clustering data from incomplete observations, which is usually an issue in multiplex imaging data that not all features are available for every data point (30). Since we chose to focus on immune complexities here, we analyzed only immune cells by using a cut-off value for CD45 (>0.07 of mean marker intensity, based on manual gating image cytometry), the pan-immune cell maker. We extracted 78,769 CD45⁺ cells from the 40 ROIs and clustered into 20 groups based on the single-cell mean intensity of 10 markers ( Figure 6 ). For our analysis, we simply explored a different number of groups, i.e., the number of subspace k (see Materials and Methods) to identify a smaller but distinct source of variation in the data with biological interpretation based on the elbow method. In general, the clustering results were comparable to the manual gating results and identified the key cell subsets including CD8⁺ T cells, FOXP3⁺ Treg cells, other T cell subsets, B cells, monocytes, and CD163⁺ myelomonocytic cells. This confirmed that our unsupervised analysis approach via SSC corroborated and complemented the manual gating strategy, which allowed the discovery of novel cell types and robust cell identification with high efficiency and objectivity from a large-scale singe-cell dataset.

We then looked at the cell density of the SSC-clustered cell groups in the disease progression ( Figure 7 , S8 ). It was noteworthy that the immune cell populations with significant change in disease progression, Treg cells and the CD163⁺ myelomonocytic cell subset, were identified using the SSC approach as group #15 ( Figures 6 , 7A ) and #12 ( Figures 6 , 7B ), respectively, and showed the similar change. Interestingly, the SSC approach also identified a myeloid cell lineage; group #7 characterized by CSF1R and CD1C ( Figure 6 ), which was not a predefined population in the manual gating strategy, also showed significant increase with disease progression ( Figure 7C ).

We reasoned that the SSC approach clustered some immune populations into distinct subsets that were not pre-defined by image cytometry manual gating. For example, groups #11, #12, and #19 were all attributed to CD163⁺ myelomonocytic cells based on the high staining density of CD163, but only group #12 was found to have a significant increase that showed the same trend in manually gated CD163⁺ myelomonocytic cells ( Figure 4C ). In addition, the cell densities of the SCC-clustered CD163⁺ myelomonocytic cell subset (Group #12) in NSQ, NDBE, Dys, and EAC ( Figure 7B ) were all proportional to the manually gated CD163⁺ myelomonocytic cells ( Figure 4C ). We therefore cautiously concluded that this specific myelomonocytic cell subset (Group #12) might be a subset of total CD163⁺ myelomonocytic cells, which is more likely to play a role in the progression of the disease compared with the other subsets. We postulated the same for the new CD1C⁺ subset (CSF1R⁺CD1C⁺, group #7) belonging to a wider CD1C⁺ population in manual gating. There were four subsets: Groups #1, #4, #8, and #10, which all showed relatively low intensities of most markers. They were labeled as unclassified immune cells and are elaborated in Discussion.

Both manual gating and SSC approaches revealed a clear increasing density of Treg cells, the CD163⁺ myelomonocytic cells, and the new CD1C⁺ cell subset. We then further studied their spatial relationship for possible cell-cell interactions during disease progression by interrogating the SSC data using two spatial analyses. First, we calculated the average shortest distance between the centers of two cell nuclei of two given cell subsets and applied 30 and 50 μm as two thresholds to evaluate cell proximity ( Figures 8A–F , S9A ). Although reflecting the absolute distance between cells, the average shortest distance could be affected by the overall cell compactness and distribution pattern. We therefore applied a second approach of neighborhood enrichment analysis that computed the likelihood that a given cell subset was neighbored by other cell subsets compared to the random background ( Figures 8G , S9B ).

We first focused on the CD1C⁺ subset and examined their neighbors. The shortest average distance observed from the CD1C⁺ subset to Treg cells followed a clear gradual decrease from NDBE to Dys and then EAC ( Figure 8A ). Especially in EAC, the majority of the cells of the two subsets were within 30 μm. Neighborhood enrichment analysis also revealed that Treg and CD1C⁺ subsets were slightly depleted in each other’s neighborhood in NDBE but significantly enriched in Dys and EAC ( Figure 8G , highlighted by yellow solid- and dotted-line boxes). A similar finding was observed in the average distance from the CD1C⁺ subset to CD8⁺ T cells, but it was generally further (>50 μm) ( Figure 8B ) compared with the distance from CD1C⁺ to Treg cells through all disease stages ( Figure 8A ). In the neighborhood of CD8⁺ T cells, there were increasingly enriched CD1C⁺ cells from NDBE to Dys and EAC ( Figure 8G , highlighted by purple dashed-line boxes); but in the neighborhood of CD1C⁺ subsets, the enrichment of CD8⁺ T cells varied between the disease stages ( Figure 8G , highlighted by purple solid-line boxes). The CD1C⁺ subset was unlikely to have direct cell–cell interaction with other T cells in all disease stages because they were not enriched in each other’s neighborhoods and the average distances were far ( Figures 8C, G , highlighted by red solid- and dotted-line boxes).

We then focused on the CD163⁺ myelomonocytic cell subset. The shortest average distance (<30 μm) was observed from CD163⁺ myelomonocytic cells to Treg in Dys ( Figure 8D ), which also had the most significant neighborhood enrichment ( Figure 8G , highlighted by green solid-line boxes in Dys). The shift of their neighborhood phenotype from a slight depletion in NDBE to enrichment in Dys ( Figure 8G , highlighted by green solid-line boxes, NDBE vs. Dys), mirrored the shift of the proportion of CD163⁺ myelomonocytic cells, which increased from 15% to 45% ( Figure 5B , NDBE vs Dys). This was consistent with literature reporting that Treg cells induce the polarization of CD163⁺ myelomonocytic cells (referred to as M2-like macrophages in the literatures) (31–33). In EAC, the average distance from CD163⁺ myelomonocytic cells to Treg cells was close ( Figure 8D , EAC), but it was unlikely that they were spatially correlated, as revealed by the neighborhood enrichment analysis ( Figure 8G , highlighted by green solid- and dotted-line boxes, EAC). This was likely due to Treg cells having significantly higher density in EAC ( Figures 4B and 7A ), such that they were spatially closer to all their neighbors but not specifically to the CD163⁺ myelomonocytic cell subset. The average distance from CD163⁺ myelomonocytic cells to CD8⁺ T cells was greater overall ( Figure 8E ) than to Treg cells ( Figure 8D ). CD8⁺ T cells were increasingly enriched in the neighborhood of CD163⁺ myelomonocytic cells in disease progression ( Figure 8G , highlighted by blue solid-line boxes), but variation was observed in the vice versa ( Figure 8G , highlighted by blue dotted-line boxes). The spatial relationship between CD163⁺ myelomonocytic cells and other T-cell subsets was unclear, where the average distance decreased with disease progression ( Figure 8F ), but the neighborhood phenotype varied ( Figure 8G , highlighted by black solid- and dotted-line boxes). It was not unexpected as this T-cell group was probably a diverse population consisting of multiple T helper cell subsets.