Resampling approaches for the quantitative analysis of spatially distributed cells

Several approaches are employed in image segmentation, ranging from traditional techniques based on the threshold values and cellular features (including size, shape and edge) to advanced methods utilizing machine learning and deep neural networks. Cellular segmentation algorithms are critical components of image analysis processes and typically involve a series of computational steps that aim to identify the background and separate cells from each other, ultimately producing precise boundaries for each cell. Indeed, these algorithms are designed to automatically identify and delineate individual cells within an image. This task can be highly challenging due to variations in cell shapes, sizes, and staining intensities. In image analysis, older software like Imagescope traditionally allows us to evaluate signal intensity or pixel count thresholds by categorizing them as weak, moderate, or strong positivity, but cannot perform cell segmentation. Over recent years, there has been a growing emphasis on cellular segmentation, driven by the advent of new methods and a heightened focus on automation and precision. Contemporary, cutting-edge software tools have evolved from foundational image processing algorithms developed decades ago, including ImageJ, Cellsegm, and CellProfiler [1,2]. HALO is a notable addition to these recent segmentation tools, which facilitates automated tissue segmentation, serial section analysis, and TMA processing (as referenced in HALO). Using HALO, we investigated how cells are spatially distributed. Moreover, we identified different cell populations by recording signal intensities.

Classifying cells in populations based on the protein markers is a common procedure in clinical medicine and scientific research. The study of those populations helps us in tissue characterization and understanding the cell interaction behaviour. A cell population could be distributed in clusters and could preferentially aggregate or segregate with other specific populations. The recorded cell set X = (x₁, x₂, …, x_N} can be described as a point pattern process within the domain De R². In this context, the D domain is affected by the presence of multiple irregularly shaped interstitial cell-devoid spaces, which correspond to anatomical features (e.g. vessels, lipid vacuoles, glandular lumina) or tissue artefacts (e.g. tissue fractures), and whose coordinates are unknown (Fig. 1). These interstitial spaces limit analyses conducted using the point pattern approach; specifically, the calculation of the domain area is not accurate. Despite some authors having utilized classical clustering measurements in spatial cytometry [3, 4, 5], no bias correction method has been introduced in the literature to account for interstitial cell-devoid spaces. To address this gap, we propose novel resampling methods for testing spatial randomness and evaluating relationships among different cell populations. Our approach doesn't need the domain area estimation and provides non-biased clustering measurements.

HALO [6] is one of the most widely adopted software tools for tissue image analysis in multiplexed staining approaches. It allows accurate quantitative tissue analyses in digital pathology, also integrating artificial intelligence [7, 8, 9]. In detail, it enables quantifying various histopathological features in digital slides with high reproducibility. The HALO platform (v3.2.1851.229, Indica Labs) image analysis tool encompasses different modules for brightfield or fluorescence microscopy. Besides the automated tissue segmentation, the different modules evaluate the spatial distribution, maintaining an interactive link between cell data and images. In our spatial cytometry dataset, we adopted a specific immunohistochemistry module for recognizing each cell's nucleus, cytoplasm, and membrane. The adopted module reports staining intensity metrics and cell coordinates within the tissue (see supplementary data).

Spatial Analysis offers the ability to identify cells’ proximity and spatial distribution in individual tissues. It determines the average distance and the number of unique neighbours between cell populations.

In this section, we present three spatial analysis methods to study spatial cell distributions and relations;

• Spatial randomness test evaluates the presence of cell clustering behavior.

• Spatial dependence test evaluates the presence of aggregation and segregation effects between two cell populations.

• Distribution equality test evaluates if two cell populations are differentially distributed over the domain.

All these methods are based on resampling approaches conceptually close to a bootstrapping procedure without replacement. Therefore, the resampling under the null hypothesis will provide the null distribution of the statistics of interest.

Let's consider a set S of N cells and a subset A of n_a. lls identified by a specific biomarker. Our approach evaluates the clustering degree of subset A by comparing the sample average nearest neighbour (NN) distance x̄_a. within the subset A with the distribution of the average NN distances under the hypothesis of spatial randomness.

The average NN distance of the subset A is equal to:

where, d_i^a. is the distance between the i^th cell of set A and its nearest neighbour cell inside A. The empirical distribution of x̄_a under the hypothesis of spatial randomness is obtained by using the following resampling procedure:

• K groups of n_a cells are randomly selected from the whole set S.

• The average NN distance of the k^th random group is calculated as.

where, d_i^k the distance between the i^th cell of set k and its nearest neighbour cell inside k.

We obtain a vector of K values that compose the empirical null distribution of x̄_a. Therefore, the p-value for spatial randomness testing can be calculated as the proportion of values from the empirical distribution lower or equal to the observed average NN distance x̄_a:

As an additional advantage of our approach, we can calculate an empirical version of the Clark-Evans aggregation index [10] as the ratio between the observed average NN distance and the average value of the null distribution:

The empirical R index has the same interpretation of the Clark-Evans index. Indeed, R < 1 indicates aggregation (i.e., clustering), whereas R > 1 indicates segregation. Finally, a value of R ≈ 1 indicates spatial randomness.

We can extend the criterion used in the previous testing procedure to evaluate if a cell population is significantly close or far from another one. Let's consider a set S of N cells and two subsets, say A and B, identified through two different biomarkers, with A ∩ B = Ø. We want to evaluate if A cells are significantly close or far from B cells. To reach this goal, we use a resampling approach that compares the average NN distance x̄_ab from A cells to B cells with the distribution of x̄_ab under the hypothesis of spatial random distribution of the B cells. The average NN distance x̄_ab is calculated as:

where, d_i^ab is the distance between the i^th cell of set A and its nearest neighbour cell of B. The empirical distribution of x̄_a under the hypothesis of spatial randomness of B cells is obtained using the following resampling procedure;

• After excluding the A cells from the S set, K groups of n_b cells are randomly selected from the remaining N-n_a cells.

• The average NN distance of the k^th random group is calculated as.

where, d_i^ab is the distance between the i^th cell of set A and its nearest neighbour cell inside k.

This resampling procedure gives a vector of K values that composes the empirical null distribution of x̄_ab.

Considering the null distribution, the dependence test can be defined as a two-tailed test or as a one-tailed test to distinguish aggregation and segregation effects. If we want to evaluate if A cells are significantly close to B cells (aggregation effect), the p-value can be calculated as the proportion of values from the empirical distribution lower or equal to the observed average NN distance x̄_ab:

On the other hand, if we want to evaluate if A cells are significantly far from B cells (segregation effect), the p-value can be calculated as the proportion of values from the empirical distribution greater or equal to the observed average NN distance x̄_ab:

Finally, the two-tailed p-value can be calculated as:

where M is the median value of the empirical null distribution. Additionally, we can calculate an interaction index through the following ratio:

where R_ab < 1 indicates an aggregation effect of A cells on B cells, R_ab > 1 indicates a segregation effect of A cells on B cells, and R_ab ≈ 1 indicates the absence of a spatial association between A and B.

It is important to underline that the spatial relation between A and B cells could be asymmetrical. Indeed, the spatial behaviour of NN cells of A could be very different from that of the NN cells of B. Therefore, the test results have to be interpreted from the point of view of the A population in terms of NN cells. In other words, we evaluate if the NN cells of A involve a significant number of B cells.

Until now, we have assumed that A ∩ B = Ø. If there are n_ab double-positive cells that belong to both A and B populations, some observed d^ab distances are equal to zero. In this case, we fix the number of double-positive cells equal to n_ab in each random sample. This strategy allows us to run the test despite double positive cells, making the testing procedure more conservative. Moreover, we suggest considering a part of the double-positive cells during the testing to understand their spatial relation with the other populations.

Let's consider two cell populations, say A and B, identified through two different biomarkers, with A ∩ B = Ø.

We propose a statistical test to evaluate the hypothesis of spatial distribution equality between A and B populations. A rejection of this hypothesis indicates a different spatial distribution of the two cell groups.

Our statistical approach is conceptually close to the quadrat count method [11]; the domain of interest is divided into m subareas (empty subareas are eventually excluded). Then, we calculate the difference between the number of the two types of cells present in each subarea. Statistical test C compares the observed differences with the expected differences under the hypothesis of spatial cell mixing:

where, d_i is the difference between the number of A cells and B cells in the i^th subarea, and $E$[d_i] is the expected difference calculated as:

where n_i is the total number of cells in the i^th subarea, and p is the proportion of A cells in the whole domain, that is, p = n_a/(n_a + n_b).

Although the previous formula looks familiar, the C variable doesn't follow a Chi-square distribution because the expected difference in the denominator is expressed as an absolute value. For this reason, we empirically calculate the null distribution through a resampling approach:

• The cell groups A and B are randomly mixed K times, fixing their sizes n_a and n_b.

• The C statistic is calculated for each randomly mixed scenario.

Therefore, the p-value can be calculated as the proportion of values from the empirical null distribution greater than or equal to the observed statistical test value c:

Like any quadrat count approach, our method is influenced by the quadrat size [12]. On the other hand, an advantage of our approach is that quadrats that go out of the domain don't bias the C statistics calculation because the expected difference d_i is calculated based on the cell number n_i in the i^th subarea.

Adopting the same criteria of the previous indexes, we can calculate a dissimilarity index as:

where, I_C ≠ 1 indicates that A and B are differentially distributed over the domain, and I_C ≈ 1 indicates that A and B cells are equally distributed.

We carried out simulation analyses to evaluate the p-value accuracy of the introduced statistical tests. Specifically, for each testing procedure, the null distribution of the test statistic has been obtained by calculating the statistic values among K two-dimensional random fields containing n_a = p · N cells of set A and N = 10.000 total cells. Then, we obtained the p-value distribution by performing K₂ additional simulations, each providing a statistical test value x_i, with i ∈ {1, …, K₂}. Therefore, the p-value p_i has been calculated by comparing the statistic value x_i with the null distribution previously obtained.

Table 1, 2, and 3 show the p-value distribution properties varying the p = n_a/N parameter. For K≥1.000, we observed a uniform distribution of p-values in almost all simulations, with a p-value accuracy close to 0.001. The results highlight high p-value reliability for all the introduced statistical tests.

As proof of concept, we applied our statistical tests to real spatial cytometry data from the area segmentation of lung and tonsil samples. In particular, we focused on lymph node germinal center (GC) compartments. The

GC is a specialized microenvironment formed within the B-cell follicles in secondary lymphoid tissues and is subdivided into two distinct functional compartments: the Dark Zone (DZ), the site of B-cell proliferation and somatic hypermutation (SHM); and the Light Zone (LZ), where antigens are presented to follicular T helper (Tfh) cells for selection [13].

Table 1 reports the spatial randomness test results from one lung sample (Fig. 3A) and two tonsil samples (Fig. 3B and 3C). Those samples clearly show a clustering behaviour of cell populations, while the bottom part of Table 1 reports the results from a randomized sample (Fig. 3D). As we expected, clustering p-values calculated on the random sample are not significant. Table 2 reports the results from the spatial dependence test. In particular, we have separately evaluated the aggregation effect of cells expressing the dark-zone associated B-cell marker AID on cells expressing the DNA-damage marker gH2AX and the segregation effect of AID+ B cells with CD3-expressing T cells. While the bottom part of Table 2 reports the results from a randomized sample (Fig. 3D). As we expected, the randomization eliminates the spatial dependence among cell populations.

Fig. 4, 5, 6, and 7 report the cell spatial distributions on which we applied the spatial equality test. In the right panel of each figure, we report the d¡ observed differences and the $E$[d_i] expected differences obtained from the quadrat count approach. In particular, Fig. 4 shows a clear cell spatial separation (p-value < 2e-04), Fig. 5 reports two cell sets artificially randomly mixed (p-value = 0.96), Fig. 6 shows two cell populations naturally mixed (p-value = 0.35), and Fig. 7 shows two cell population spatially separated (p-value = 0.0016).

Our results provide a quantitative layout for interpreting specific biological features of the investigated tissue. Specifically, we demonstrated a spatial association between AID-expressing B cells of the dark zone, which undergo replication and mutational activity, with the DNA-damage-associated marker pgH2AX. Moreover, we could support by quantitative analysis, consistent segregation of the AID+ DZ cells with T-cells residing in the germinal center, a feature that prompts intriguing hypotheses regarding the selective regulation of T-cell presence in the DZ.

This paper introduces a novel testing procedure for spatial cytometry data. The procedure is not affected by the presence of interstitial cell-devoid spaces in the domain, which typically determines a bias in the outcome of classical clustering measurements.

Even though certain authors have applied traditional point pattern techniques in spatial cytometry [3], no method has been proposed in the literature to consider the presence of interstitial cell-devoid spaces. Those spaces don't allow a reliable domain area estimation, causing inaccuracy in clustering measurements and p-value calculation. Our research significantly advances the existing literature by providing a point pattern approach that avoids the domain area estimation.

The presented methods can be used to test spatial randomness, aggregation effects, segregation effects, and spatial distribution equality, highlighting the spatial relationships among cell populations. Spatial interaction and clustering indices can be estimated to describe spatial cell behaviour.

Our methods can be applied to any point pattern analysis in which an inaccurate estimation of the domain area limits the study of sub-populations. Indeed, the elements that don't belong to the group of interest can be considered to perform the resampling and construct the null distributions.

The simulations support a high p-value accuracy for resampling numbers of 1,000 or more. The results obtained for real data provide a clear view of spatial cell behaviour, while the analysis of artificially randomized data yields not-significant p-values coherently with the null hypothesis.

This article's Supplementary data and R functions can be found on the SpaceR GitHub page (https://github.com/GBertolazzi/SpaceR). The SpaceR software allows the reader to use the R functions for analyzing his own data. The supplementary dataset can be used as an example to run the analysis reported in the manuscript.

Claudio Tripodo formulated the research problems. Giorgio Bertolazzi developed and implemented the method under the supervision of Michele Tumminello and Claudio Tripodo. Gaia Morello and Beatrice Belmonte collected the data. Giorgio Bertolazzi and Michele Tumminello analyzed the data. All the authors wrote and revised the manuscript.

This study has been supported by the Italian Ministry of Education, University and Research (MIUR) through the “PON Research and Innovation 2014-2020” to G.B.; by the National Biodiversity Future Center (NBFC) CN00000033 (CUP UNIPA B73C22000790001), and through the OBIND project N.086202000366 (CUP G29J18000700007) to M.T.; by the Italian Foundation for Cancer Research (AIRC) through the 5x1000 I.D. 22759 Grant and AIRC Accelerator Award ID. 24296 to C.T.; and by the Italian Ministry of Education, University and Research (MIUR) Grant 2017K7FSYB to C.T.