1 Introduction

Pancreatic cancer is related to be the second leading cause of the cancer-related death in the global world by 2030 [1]. Despite the fact that there have been continues advances, the most prevalent type of pancreatic cancer, Pancreatic ductal adenocarcinoma (PDAC), is still one of the most devastating tumor, with a poor 5-year survival rate of approximately 10% in the United States [2, 3]. Immunotherapy, such as therapies targeting immune checkpoint molecules, has emerged as one of the most promising treatment options for various solid tumors [4]. However, despite the success observed in other cancers, progress in improving outcomes for PDAC patients has remained limited [5]. The part of reasons may be due to the tumor microenvironment (TME), which is the primary key to the poor clinical responses to immunotherapy. Studies pointed out that the TME not only promotes the development of tumor but also serve as a formidable barrier to the effectiveness of immunotherapies by diminishing the efficacy of the drug [4, 6,7,8].

Recent advances in high-dimensional multiplexed imaging techniques have substantially deepened our understanding of the TME in PDAC. The TME in PDAC comprises cancer cells, stromal cells, and immune cells. However, the contribution of immune cells to tumor progression and therapeutic resistance warrants careful examination. Previous studies highlight the significance of the spatial landscapes of various immune populations, including B cells, T cells, myelomonocytes, and cancer-associated fibroblasts, in relation to prognostic outcomes [9,10,11,12,13]. Despite this, the spatial location of immune cells can significantly impact patient survival [14,15,16]. Therefore, it is crucial to conduct separate analyses based on the spatial localization of cells. Among these immune cells, T cells have been explored for their role in antitumor responses, however, despite their recognized presence, the specific role of B cells remains less well understood [17].

The gap in understanding has sparked growing interest in investigating the role of B cells across various tumor types. Until now, many researchers have started to pay more attention to the functions of B cells of different types of tumors. For triple-negative breast cancer, spatially distributed B cells indicate good prognostic tumors, revealed by Wortman et al. [6]. Some studies reported the controversial results of the effect of B cells in gastric cancers [18, 19]. For instance, Huifang Liu et al. found the considerable association between dense peritumoral TIGIT+CD20+ B cells infiltration and poor clinical outcomes [20]. In contrast, observed by Audrey et al., a great number of CD20+ B cell follicles bring patients better survive [21]. There have been limited investigations into the relationship between B cells and survival in pancreatic carcinoma patients. Jiaren Xu and colleagues demonstrated that B cells exhibit markedly different functions within tertiary lymphoid structures compared to the matrix microenvironment [22]. Some studies suggested that elevated B cell levels and increased B-cell aggregates are associated with improved prognosis [23,24,25]. Conversely, other researchers highlighted that the prominent presence of B cells contributed to worse outcomes by promoting pancreatic tumor growth [26,27,28]. Beyond their involvement in forming tertiary lymphoid structures (TLS), B cells dispersed throughout the TME also fulfill critical roles, including antibody production, cytokine secretion, and antigen presentation to helper T cells. Given this gap, we suppose the controversial results above are due to the spatial heterogeneity of B cells. To be specific, the spatial pattern of B cells within the TME differentially impacts tumor growth and patient survival.

Traditional methods, such as light microscopes, are not able to provide precise, quantitative, and replicable research findings. Nevertheless, supported by digital pathology, which is an increasingly popular subject, pathologists can tract information of interests from the digital pathology slides [29, 30], which offers a comprehensive and enriched view. In this study, we designed and implemented an automated workflow which utilized several key metrics to capture spatial heterogeneity. We divided the tumor into three distinct regions. In each of these regions, spatial distribution analysis and a density-based clustering algorithm were applied to understand the spatial pattern of CD20+ B cells in relation to patient survival. Through this comprehensive and computer-based workflow, we effectively reduced the need for manual intervention, thereby increasing the efficiency and reproducibility of our analysis.

2 Methods and materials

The input data for our analysis consists of digital pathology slides stained via immunohistochemistry (IHC) to label CD20+ B cells in post-surgical PDAC patients who did not receive any adjuvant therapy. For each patient slide (n = 45), we employed image segmentation to extract the coordinates of CD20+ B cells. The spatial information of these cells was analyzed using two primary methods: spatial analysis and clustering analysis. In the spatial analysis, the tumor was divided into three distinct regions, such as invasive front (IF), cancer center (CT), and cancer island (CI). There are various spatial metrics, such as occupancy rate (as area under curve, occupancy AUC), fractal dimension (FD), which were computed for each region. Meanwhile, in cluster analysis, B cell clusters were identified using the HDBSCAN algorithm and metrics such as cluster density, and coverage ratio were evaluated. These analyses allowed us to quantitatively assess the spatial distribution of B cells within the TME. The results from both spatial and cluster analyses were correlated with patient survival outcomes to evaluate the prognostic relevance of B cell spatial pattern. The overall workflow is illustrated in Fig. 1.

Fig. 1
figure 1

Overall workflow of spatial pattern analysis of CD20+ B cells Overall workflow starts with the preparation of whole slide images (WSI) and pathologist annotations. Tumor of each slide was defined as specific tumor regions, including IF, CT and CI. Following this, image segmentation is performed to extract the coordinates of CD20+ B cells from each slide. These cell coordinates are then projected onto the defined tumor regions to construct a comprehensive cell map for subsequent spatial analysis. The spatial pattern analysis includes two main components: spatial distribution analysis and density-based clustering analysis. In spatial distribution analysis, the occupancy rate and fractal dimension are used to assess the spatial characteristics of CD20+ B cells. In the clustering analysis, B cells are grouped into clusters, followed by the calculation of key metrics including cluster density and coverage ratio. Finally, all derived spatial metrics within TME are subjected to survival analysis to assess their prognostic significance. This workflow is repeated across all patient samples to capture the association between spatial pattern of CD20+ B cells and patient prognosis outcomes in three tumor regions

Full size image

2.1 Study populations

Forty-five patients who underwent surgical resection for PDAC at the First Affiliated Hospital of Fujian Medical University from 2010 to 2021 were included. None of the patients had undergone neoadjuvant and radio therapy before surgery. Formalin-fixed paraffin-embedded blocks of each PDAC specimen were cut into 4-μm-thick slices and stained with H&E. All original H&E slides were reviewed, and a representative slide that included both the tumor tissue and the adjacent healthy tissue was selected. The histological diagnoses and grades were confirmed according to the latest World Health Organization criteria [31]. Patients’ demographic and clinical details, including age, pN stage, pT stage, and outcomes, were retrieved from their medical records. OS was defined as the interval from the date of the initial diagnosis to the date of death from any cause. Patients who died within 1 month after surgery were excluded.

2.2 Immunohistochemical staining

One representative 4-μm-thick slide that contained tumor and adjacent paracancerous tissue was selected from each patient by two pathologists. All the slides were stained for CD20 (C8/144B, 1:100, Dako, Denmark). Whole slides were prepared for IHC using monoclonal antibodies according to the manufacturer’s protocol. After de-waxing, hydration, and antigen retrieval (EDTA at pH 8.0, 37 °C for 30 min) and endogenous peroxidase blocking (3% H2O2 at room temperature for 10 min), a Dako EnVision FLEX + detection system (DK-2600, Dako, Denmark) and an Ultra view Universal DAB Detection Kit (Ventana, AZ, US) were performed. Human tonsil tissues were used as the positive control, and phosphate buffered saline (instead of the primary antibody) was used as the negative control.

2.3 Image process

To perform spatial analysis, we extracted the coordinates of CD20+ B cells from whole slide images (WSI). Due to physical memory constraints, each WSI was down sampled by a factor of four and subsequently divided into patches with a resolution of 400 × 400 pixels using Groovy scripts. Following this, QuPath (v3.5.0) which is an open-source software known for its robust technical support, effective version control and powerful capabilities in digital pathology image analysis, was used to segment CD20+ B cells within each patch [3, 32]. A wide range of published studies have implemented an trained unsupervised watershed algorithm, which was used to detect positive(stained) nucleus in our analysis [3, 33,34,35,36,37,38,39]. During this process, a high threshold was applied to remove background noise and nonspecific immunostaining. Once detection has finished, the coordinates of the positive cells within each patch were recorded, extracted and saved as corresponding txt files for further analysis. To evaluate the performance of detection, we conducted a thorough validation analysis using 20 WSI. Samples (N = 30) were taken from each slide using random sample method in the three regions: IF, CT and normal tissue (N). This approach resulted in a total of 600 subregions for validation, providing a robust and diverse dataset to assess classifier performance across heterogeneous tissue environments. To evaluate the QuPath classifier, we performed both manual annotations (by an experienced pathologist) and automated detection. Cells detected by both manual annotation and software are classified as true positives (TP). In contrast, cells identified by the software but not confirmed by manual annotation are classified as false positives (FP). Additionally, cells missed by the software but identified manually are considered false negatives (FN). The Sensitivity (recall rate, R) and Precision (P) of the software algorithm and the standard error are presented as:

$$ R = \frac{TP}{{TP + FN}} $$
(1)
$$ P = \frac{TP}{{TP + FP}} $$
(2)
$$ se_{R} = \sqrt {\frac{R(1 - R)}{{TP + FN}}} $$
(3)
$$ se_{P} = \sqrt {\frac{P(1 - P)}{{TP + FP}}} $$
(4)

2.4 Region characterization and Cell projection

For each case, we read the coordinates of the detected positive cells within each patch from the corresponding files. These scattered cells represent the local coordinate system. After that, we transformed the local coordinate system into a global coordinate system and generated a two-dimensional map where the coordinates were subsequently projected. Simultaneously, the whole tissue outline and the boundary between cancer tissue and normal tissue were selected and annotated by pathologists in ASAP (v2.2), after which they were converted into coordinate sets. The tissue outline is closed so that the coordinate set of it can form a polygon, which is the whole tissue. The boundary was buffered 0.5 mm inwards and 0.5 mm outwards to create IF region, which defines the tumor boundary and divides the whole tissue into CT and normal tissue [3, 35, 40,41,42]. We further defined the cancer island (CI) as consisting of both CT and IF.

2.5 Quantitative Metrics for Spatial Analysis

To analyze the spatial pattern of cells, we quantified the B cells within IF, CT, and CI regions respectively. We utilized python’s ‘shapely’ library to generate minimum bounding rectangles for each region. Afterward, the IF, CT, and CI areas were fully covered using square grids of varying sizes, which made sure that three regions were comprehensive partitioned into subregions for the application of the following statistical metrics. The grid size was chosen to be small enough to capture variations in cell density and large enough to examine concentrated spatial distribution pattern. In this case, the minimum grid size was set to 10 µm (\(L\) = 10 µm), and the maximum size was set to 600 µm (\(L\) = 600 µm).

2.5.1 Occupancy analysis

Occupancy analysis allows for the assessment of the spatial distribution of positive cells. For each region, a square gird of side length \(L\) is used, with the area of each square being \({L}^{2}\). The total number of square grids is denoted as \(n({L}^{2})\). In each square grid (subregion), a detection process was carried out to determine the presence of B cells. Based on the results, if B cells were detected within the gird, it was marked as ‘1’; otherwise, it was marked as ‘0’. The occupancy rate (OR) was defined as the proportion of girds marked as ‘1’ relative to the total number of subregions, represented by the following formula:

$$ OR = \frac{{n(L_{M = 1} )}}{{n(L_{M} = 1) + n(L_{M = 0} )}} $$
(5)

As the side length \(L\) of the grid increases, the occupancy rate is calculated to describe the spatial distribution of B cells at different scales of \(L\). It is important to note that occupancy is influenced by the mean density of cells, as a higher average cell density increases the probability that a grid will be labeled as ‘1’.

2.5.2 Fractal dimension

The fractal dimension provides a measure of spatial occupancy and complexity within a given region. There are several methods for estimating fractal dimensions, including box-counting, mass-radius, dividers, and Bouligand-Minkowski techniques. In this study, we employed a variant of the box-counting method to represent the spatial distribution of cells using the number of square subregions marked as ‘1’, denoted as \(n({L}_{M=1})\) [43]. This approach facilitates the straightforward estimation of the fractal dimension of tissue structures. For each grid with a side length \(L\), the area of the window is \({L}^{2}\), and the number of square windows \(n({L}^{2})\) is inversely proportional to \(1/{L}^{2}\). Similarly, the number of square subregions marked with ‘1’, \(n({L}_{M=1})\), is proportional to \(1/{L}^{d}\), where \(d\) \(\le\) 2 [44,45,46,47,48,49]. To quantify the analysis dimension across a range of scales, we defined two ranges: a small scale (10–40 µm) and a large scale (200–600 µm). If the tissue exhibits self-similarity at a particular scale, it should adhere to the following power law: \(n({L}^{2})=A/({L}^{d})\), where \(A\) is a constant influenced by \(n({L}^{2})\). To minimize the impact of tissue size on \(n({L}^{2})\), we defined:

$$ FD = S(L^{2} ) = \frac{{d[\log \;n(L^{2} )]}}{d[\log \;(A/L)]} $$
(6)

FD is calculated as the slope S of the line fitted to the points using a least-squares linear regression.

2.5.3 Difference in Fractal Dimension

The difference in fractal dimensions (FD difference, ΔFD) across different scales is used to quantify the spatial clustering of cells within the tissue. This measure is defined as:

$$ \Delta FD = FD\;(l\arg e\,scale) - FD\;(small\;scale) $$
(7)

where FD (large scale) represents the FD at a larger scale, and FD (small scale) represents the FD at a smaller scale [50]. It quantifies the shift in the spatial distribution pattern of cells across these scales, offering insight into how B cells are spatially organized within the tumor microenvironment [51]. A large ΔFD indicates that at larger spatial scales, B cells exhibit a dispersed, two-dimensional plane-like distribution, while at smaller scales, they demonstrate a more compact, clustered arrangement. This shift reflects a transition from localized clustering to a more spatially spread-out distribution as the observation scale increases. Conversely, a small ΔFD suggests that the spatial distribution of B cells is consistent across both large and small scales, indicating a stable, self-similar clustering pattern. In this scenario, B cells maintain a cohesive, localized aggregation at multiple scales, signifying a persistent clustering behavior.

2.5.4 Cluster analysis

To investigate the impact of cell clusters on patient prognosis, we employed the HDBSCAN algorithm from the python’s ‘hdbscan’ library to cluster B cells and identify immune cell clusters. HDBSCAN is a density-based clustering algorithm that performs global point pattern analysis. The algorithm operates through a series of steps: building the minimum spanning tree, constructing the cluster hierarchy, condensing the cluster tree, and extracting the clusters [52, 53]. For the parameter settings, we set the ‘min_cluster_size’ to 30, which determines the smallest size of clusters to be considered. Additionally, we set the parameter ‘min_samples’ (clustering stringency level) to 15. We measured each cluster area as the area of concave hull, and decided to define two statistic indicators, which are cluster density and coverage ratio:

Cluster density: Measures the number of B cell clusters within a given region. Cluster density is defined as:

$$ D\;(cluster) = \frac{n\;(cluster)}{{A\;(region)}} $$
(8)

where \(n(cluster)\) represents the number of clusters and \(A(region)\) denotes the area of the region.

Coverage ratio: Measures the proportion of the region occupied by B cell clusters. Coverage ratio is defined as:

$$ R\;(cluster = \frac{\sum A (cluster)}{{A\;(region)}} $$
(9)

where \(\sum A(cluster)\) is the sum of the areas of each cluster within region and \(A(region)\) is the area of the entire region.

2.6 Statistical analysis

All statistical analyses were performed using python. Patients were grouped based on the mean value of spatial metrics and the survival analysis was performed using the Kaplan–Meier method combined with the log-rank test. Mann–Whitney U Test was used to evaluate the significant difference between LTS and STS within statistical metrics. The area under curve (AUC) was used to facilitate comparison of the occupancy rate. A p value of < 0.05 was considered statistically significant.

3 Results

3.1 Characteristics of the study populations

The study included 45 PDAC cases, comprising 27 male and 18 female patients with a mean age of 59.6 years (28–86 years). The median tumor size was 3.96 cm (range, 0.7–9.0 cm). Regarding the tumor histological grade, there were 13 cases of well-differentiated tumors, 9 cases of moderately differentiated tumors, and 23 cases of poorly differentiated tumors. For TNM staging, 7, 23, and 15 cases were classified as pT1, pT2, and pT3; 22, 18, and 5 cases were classified as pN0, pN1, and pN2; and 43 and 2 cases were classified as pM0 and pM1, respectively. There were 34 patients at TNM stage I–II and 11 patients at TNM stage III–IV. Follow-up data were available for all patients, among whom 32 died of PDAC. The avarage OS of all patients was 16.69 months (range, 2–84 months). To perform a differential analysis based on OS, the patients were divided into two groups: long-term survivors (LTS), defined as those with OS longer than 24 months, and short-term survivors (STS), defined as those with OS shorter than 24 months. There are 17 patients in the LTS group, while the STS group included 28 patients (Table S1). The advantage of this splitting strategy is that it could maximize the differences in spatial heterogeneity and identify the factors that affect prognosis outcomes [54]. The overall survival of 45 patients was presented in Table S2.

3.2 Segmentation results

After segmentation, we assessed performance by comparing manual identification with automated detection across 30 subregions from 20 WSI. The classification results were as follows: TP = 7002, FP = 1625, FN = 559. With these values, we calculated the performance metrics: Recall rate \(R(\%)=92.6\pm 0.3\), Precision \(P(\%)=81.2\pm 0.4\). To further validate the consistency of the QuPath classifier with manual annotations, we calculated the Spearman’s rank correlation coefficient between the classifier predictions and manual cell counts. The resulting correlation coefficient was 0.984, indicating a very strong positive correlation. The classifier’s performance across the three distinct tissue regions (IF, CT, and N) is summarized in Table 1, which details the region-specific metrics to illustrate the robustness of the QuPath classifier in varying tissue environments. Additionally, Fig. 2 provides a visual comparison between the classifier’s segmentation and manual annotations. To further demonstrate the generalizability of the classifier, we randomly selected 6 subregions from the IF, CT, and N regions to assess its performance.

Table 1 Classifier’s results across invasive front, cancer center and normal tissue
Full size table
Fig. 2
figure 2

Visual Comparison between the Classifier’s Segmentation (red outlines) and Manual Identification (green dots). 600 subregions are selected from 3 different regions across slides to evaluate the performance. a, b Sample 29 and 33 from invasive front. c, d Sample 117 and 231 from cancer center. e, f Sample 79 and 171 from normal tissue

Full size image

3.3 Data process and statistical characteristics

As illustrated in Fig. 3a, each patch represents a subregion, and the WSI was segmented into multiple subregions for detection. Red outlines indicate detected positive cells (CD20+ B cells), while blue outlines represent detected negative cells. Figure 3b illustrates cell visualization, clearly distinguishing positive and negative cells using red and blue outlines, respectively. After detection, the coordinates representing each region were converted into polygons. These regions were color-coded, with bright yellow representing IF, peach orange for CT, and the combined IF and CT regions marked as CI with a red outline (Fig. 3c, d). Following this, the cell coordinate system was reconstructed, as shown in Fig. 3e, by mapping local cell coordinates back to their respective subregions. As demonstrated in Fig. 3f, we projected the regions and cells together, enabling clear separation of B cells within the three distinct areas.

Fig. 3
figure 3

Cell detection, projection, and clustering Process a An original whole slide image with segmentation lines, with the segmented patches used for detection. b Cells visualization and cells detection results of the yellow-boxed subregion highlighted in panel (a). In this process, CD20+ B cells are detected with red outlines. c An original whole slide image for region characterization and cell projection. d Using region coordinate sets to form three regions, invasive front (IF, yellow region), cancer center (CT, peach orange) and cancer island (CI, red outline). e A cell map which is reconstructed based on the coordinates from detected cells. f Projection result, formed by combining three regions from panel (d) and cell map from panel (e), where CD20 + B cells are represented in dark red. g An original whole slide image for clustering. h Clustering result. Each color corresponds to a detected cell cluster. Concave hulls (colored dashed line) were calculated for further studies. The dark red points represent CD20 + B cells, while the yellow region indicates IF, the peach-orange region represents CT, and the red outline marks CI

Full size image

We observed that B cells were distributed across different regions, manifesting in both clustered and dispersed patterns. To further analyze B cell clustering, we applied the HDBSCAN algorithm, resulting in the clustering output shown in Fig. 3g, h. The concave hulls of different colors represent distinct clusters, with varying degrees of cluster density. For a more detailed regional analysis, we classified clusters based on whether their centroids were located within a specific region. For clusters spanning multiple regions, we included them in the count for each relevant region.

As shown in Table 2, in terms of occupancy AUC, IF exhibited the highest AUC. Significant differences were observed between IF and both CT and CI (p < 0.001), while no significant difference was noted between CT and CI. For ΔFD, the values across IF, CT, and CI were similar, with no statistically significant differences (p > 0.05). Cluster Density was highest in IF, followed by CI, with the lowest density observed in CT. The differences between IF and the other regions were highly significant (p < 0.001), and the difference between CT and CI was also significant (p = 0.016539). Coverage Ratio exhibited the most variability in IF, while CT and CI exhibited lower and more consistent values. A significant difference was found between IF and CT (p = 0.000884), whereas the differences between CI and both IF and CT were not significant.

Table 2 Statistical characteristics of invasive front, cancer center, and cancer island
Full size table

3.4 Occupancy rate of CD20+ B cells was not significantly correlated with OS in PDAC

We hypothesized that the spatial pattern of CD20+ B cells could provide valuable insights into TME. To evaluate this, we conducted occupancy analysis in three distinct tumor regions. Figure 4a visualized the occupancy AUC values for each case across the three regions. Generally, the AUC values increased simultaneously across IF, CT, and CI. Figure 4b depicted the average occupancy rates of B cells across three regions for 45 cases, with square grids of varying lengths \(L\). As grid size increased, the average occupancy rate also increased, reflecting the larger area where more B cells were captured by each grid. The red, yellow, and blue lines represented IF, CT, and CI, respectively. By comparing LTS (solid lines) and STS (dashed lines), we observed that LTS tends to have lower occupancy rates across all regions, implying a more spatially concentrated B cell distribution compared to STS. Figure 4c–e illustrated the comparison of occupancy AUC between LTS and STS groups in IF, CT, and CI, respectively. Across all three regions, the AUC values were generally smaller for LTS, though the differences were not statistically significant. According to the results of the survival analysis (Fig. 4f–h), no significant association was found between occupancy AUC and patient survival in different regions (IF: p = 0.63; CT: p = 0.61; CI: p = 0.32). Based on these findings, it could be inferred that the occupancy rate of CD20+ B cells in IF, CT, and CI was unlikely to be associated with patient prognosis.

Fig. 4
figure 4

Occupancy Analysis of CD20 + B Cells in Different Regions a Visualization of occupancy AUC values for each case in IF, CT and CI. b Line chart showing the occupancy rate for long-term survivors (solid lines) and short-term survivors (dashed lines) in IF (red lines), CT (yellow lines) and CI (blue lines). c–e Boxplots showing no significant differences for occupancy AUC of CD20.+ B Cells for long-term survivors and short-term survivors in IF, CT and CI. f–h Overall survival analyses showing no association between occupancy AUC with patient survivals in IF, CT and CI (p > 0.05). P-values in c–e are determined by Mann–Whitney U Test. P-values in f–h are determined by log-rank test. ns: not significant; *p < 0.05

Full size image

3.5 Less spatially dispersed B cells in IF region were associated with better survival outcomes in PDAC patients

FD analysis provides a means to quantify the spatial distribution of CD20+ B cells, which represents different levels of tissue observation. The ΔFD was calculated as the difference between the FD values at large and small scales, reflecting the spatial aggregation state of CD20+ B cells. A large ΔFD indicated that CD20+ B cells were more spatially dispersed, while a small ΔFD suggested that they were more spatially confined. Figure 5a, b, and c illustrated the ΔFD values in IF, CT, and CI regions for 45 samples, respectively. The ΔFD values ranged from -1 to 1, indicating variability in the spatial distribution pattern of CD20+ B cells across different cases. As shown in Fig. 5d, compared to STS, a significant lower in ΔFD for LTS was observed in IF (p = 0.05), highlighting distinct spatial behaviors. However, this trend was not observed in CT and CI region (Fig. 5e, f). The survival analysis presented in Fig. 5g demonstrated that a lower ΔFD in IF was significantly associated with improved overall survival (p = 0.05), indicating that clustered CD20+ B cells in IF contribute positively to patient outcomes. However, this pattern did not hold for CT and CI, as shown in Fig. 5h and Fig. 5i, where ΔFD had no significant impact on patient survival.

Fig. 5
figure 5

Fractal Dimension Analysis of CD20+ B Cells in Different Regions (a-c) Visualization of ΔFD values for each case in a IF, b CT, c CI. d–f Boxplots showing the FD analysis results for long-term survivors and short-term survivors in d IF, e CT, f CI. ΔFD is significant smaller for long-term survivors vs. short-term survivors in IF (p = 0.05). g Overall survival analysis illustrating the clinical significance of ΔFD of CD20+ B cells in IF (p = 0.05), using a cutoff value of 0.23 to distinguish patients with good survivals. h Overall survival analysis illustrating the clinical significance of ΔFD of CD20+ B cells in CT (p = 0.91), using a cutoff value of 0.14 to distinguish patients with good survivals. i Overall survival analysis illustrating the clinical significance of ΔFD of CD20.+ B cells in CI (p = 0.17), using a cutoff value of 0.14 to distinguish patients with good survivals. P-values in d–f are determined by Mann–Whitney U Test. P-values in g–i are determined by log-rank test. ns: not significant; *p < 0.05

Full size image

3.6 A higher density of B cell clusters and smaller coverage ratio in IF region was associated with prolonged survival in PDAC patients

HDBSCAN was applied to each case to identify B cell clusters. Figure 6a showed the cluster density values, while Fig. 6b displayed the coverage ratio values in IF, CT, and CI for 45 samples. No significant differences in cluster density were observed between LTS and STS groups across all three regions (Fig. 6c). As shown in Fig. 6d, in contrast to STS group, a significant lower in coverage ratio for LTS was found in IF (p = 0.03), whereas no significant differences were found in CT and CI. Kaplan–Meier survival analysis exhibited patients with higher cluster density in IF had significantly favorable prognosis (Fig. 6e; p = 0.003), while cluster density in CT and CI had no significant impact on patient outcomes (Fig. 6f, g; CT: p = 0.72; CI: p = 0.43). From Fig. 6h, it was evident that a lower coverage ratio in IF was significantly associated with better patient outcomes (p = 0.02). Notably, a similar trend was observed in CT (Fig. 6i), although it did not reach statistical significance (p = 0.07). In contrast, coverage ratio in CI showed no association with patient prognosis (Fig. 6j; p = 0.41).

Fig. 6
figure 6

Cluster Analysis of CD20+ B Cells in Different Regions a Visualization of cluster density values for each case in IF, CT and CI. b Visualization of coverage ratio values in IF, CT and CI. c Boxplot displaying the cluster density of CD20+ B Cells for long-term survivors and short-term survivors in IF, CT and CI. d Boxplot showing the coverage ratio of CD20+ B Cells for long-term survivors and short-term survivors in IF, CT and CI. Coverage ratio is significant smaller for long-term survivors vs. short-term survivors in IF(p = 0.03). e–g Overall survival analyses showing the clinical significance of cluster density of CD20+ B cells in e IF, f CT, g CI. Cluster density showing significant clinical difference in IF (p = 0.003), using a threshold of 0.34 clusters/mm2. h–j Overall survival analyses showing the clinical significance of coverage ratio of CD20.+ B cells in h IF, i CT, j CI. Coverage ratio showing significant clinical difference in IF (p = 0.02), using a threshold of 0.14. P-values in c, d are determined by Mann–Whitney U Test. P-values in e–j are determined by log-rank test. ns: not significant; *p < 0.05

Full size image

4 Discussion

In this study, we employed a systematic, computer algorithm-based workflow to evaluate and analyze CD20+ B cells spatial pattern in WSI from 45 postoperative PDAC patients who did not receive any preoperative adjuvant therapy. The digital processing approach was designed to quantify the spatial distribution pattern of B cells within the TME and to explore their impact on patient survival outcomes. To ensure a meaningful comparison of the immune microenvironment, we distinguished three distinct tumor regions including IF, CT and CI regions for each sample, where an automated WSI processing method was employed with statistical metrics. This study presented a reproducible and quantifiable approach for analyzing cells’ distribution pattern within TME. It minimized manual steps and addressed the challenges of limited clinical samples. Although our research only focuses on B cells in PDAC, it would provide a spatial analysis framework for digital pathology that can be extended to studies involving multiplex immunohistochemistry and immunofluorescence across various tumor types, which enables deeper and more comprehensive analysis.

In the presented study, we focused on the relationship between CD20+ B cells within TME of PDAC and patient prognosis outcomes. Given that previous studies have indicated that B cell density alone did not significantly impact prognosis [55], we introduced spatial metrics such as occupancy rate, ΔFD, cluster density, and coverage ratio, which were applied into three tumor bed: IF, CT, and CI. By employing the area under the curve to statistically compare the occupancy rate, we observed B cells were distributed across different regions, manifesting in both clustered and dispersed patterns. Compared to features of B cell spatial distribution in CT region, B cells in IF region exhibited higher occupancy rate, cluster density and coverage ratio, which indicated that spatial characteristics of B cells were associated with the location in PDAC. In lined with previous studies, some study pointed out that the TME of pancreatic cancer is related to the location of the sample, with differences in immune composition across distinct regions [56]. In addition, our results revealed that the differences between the spatial distribution of infiltrating immune cells and the locations exhibited differentiated functions in the terms of prognosis in PDAC patients. Thus, it is pivotal to specify the location of the examined sample in future research instead of simply mentioning tumor vs normal tissue.

By correlating spatial metrics across regions with patient survival, we found that in IF regions, B cells play a pivotal role in shaping patient prognosis. We observed that the smaller ΔFD of B cell distribution in IF region was significant positive associated with longer survival, suggesting a more compact spatial dispersion of B cells be able to enhance antitumor effect. Furthermore, Kaplan–Meier survival analysis also showed that patients with higher B cell cluster density and lower coverage ratio instead of a more dispersed distribution of B cells in IF region harbored significantly longer OS in PDAC, further supporting the biological role of compact clustered B cells ongoing with a protective antitumor response. Aligning with previous findings, we also found the more scattered distribution of B cells in TME did not exhibit their antitumor function in PDAC. The observation suggests that the role of B cells in PDAC is more intricate than previous understanding, which is very different from what observed for T cells, whose localization within TME is critical to predict better prognosis [14,15,16]. In most cases, the compact clustering of B cells may be related to B cell aggregates as those that found in TLS, which are ectopic lymphoid structures to enhance interactions with T cells, facilitating immune synapse formation and signal transduction, ultimately contributing to improved prognosis for PDAC patients [10, 57]. In such situations, B cells could recognize and internalize specific antigens via their receptors, subsequently presenting them to helper T cells and promoting the activation of cytotoxic T cells [58]. Nevertheless, for scattered B cells in TME, some studies suggested that the majority of scattered B cells may belong to the certain B cell subsets [26,27,28], which may contribute to tumor progression by facilitating angiogenesis and differentiating into regulatory B cells (Bregs) by releasing immunosuppressive cytokines such as IL-10, IL-35, and TGF-β to promote an immunosuppressive environment [59,60,61].

It is important to note that there are several limitations in our study. Although the presence of non-dispersed B cells and B cells clusters in IF are predictive of better outcomes of patients in PDAC in a relatively small sample, further validation should be explored in a larger independent cohort or in a multicentric study. Additionally, our study focused solely on the impact of B cells, ignoring the roles of other immune markers and their interactions. Future research should incorporate fluorescent staining samples to broaden the scope of the study. This research underscores the tremendous potential of computer-aided workflows in pathological analysis. If feasible, integrating artificial intelligence methods could further expand the current workflow, such as introducing MedSAM [62] to enhance cell recognition performance and using graph neural networks to deeply analyze graph structures in digital pathology slides.

5 Conclusion

In summary, this study aimed to quantify the spatial pattern of CD20+ B cells within the TME of PDAC patients by developing a computer algorithm-assisted spatial information quantification analysis method, capturing its relationship with patient survival across IF, CT, and CI. We figured out that the spatial distribution of B cells in IF is significantly associated with PDAC patient survivals. The presence of non-dispersed B cells and B cells clusters in IF are clinically positive to prognostic outcomes of patients. As for results in CT and CI, although it did not definitively establish the clinical significance, it demonstrated a discernible trend. The methodology shows substantial potential for applying computer-aided techniques in pathological analysis, which could help pathologists improve prognosis and guide rehabilitation therapy more effectively. However, further validation in a larger independent cohort or in a multicentric study could be taken into account.