This study included 11 treatment-naïve patients with biliary tract adenocarcinoma, representing a spectrum of anatomical subtypes: iCCA (n = 4), pCCA (n = 5), dCCA (n = 1), and GBC (n = 1). Formalin-fixed, paraffin-embedded (FFPE) tumor samples were obtained from surgical resections prior to the initiation of systemic therapy. Molecular profiling was performed on available FFPE samples using next-generation sequencing (NGS). The FoundationOne CDx assay, covering 324 genes and providing tumor mutational burden (TMB) and microsatellite stability status, was used for patients P01, P02, P03, P05, P07, P09, P10 and P11. Patient P06 was analyzed using the Myriapod NGS Cancer Panel, covering 56 genes and providing microsatellite stability status, with orthogonal validation by fluorescence in situ hybridization (FISH). No molecular analysis was available for P04 and P08.

Clinical data were collected retrospectively. All patients received first-line systemic treatment for advanced disease, consisting of cisplatin and gemcitabine with or without durvalumab (CGD or CG, respectively), according to clinical eligibility and regimen availability at the time of treatment. Overall survival (OS) was defined as the time from initiation of first-line therapy to death from any cause, while progression-free survival (PFS) was defined as the time from treatment initiation to disease progression or death. Treatment response was assessed according to RECIST (version 1.1) criteria. Patients achieving complete response (CR) or partial response (PR) were classified as responders (R), whereas those with stable disease (SD) or progressive disease (PD) were classified as non-responders (NR).

Baseline clinical and genomic characteristics are reported in Table 1, while follow-up and outcome data are summarized in Table 2.

We tested a panel of 4 antibodies using multiplex immunofluorescence (mIF) staining using a commercially available kit (Opal 6- plex Detection Kit, catalogue number NEL811001KT, Akoya Biosciences Inc., Marlborough, MA, USA) with a validated antibody panel: CD4 (Opal 480 MSI), CD8 (Opal 520 Cy5 MSI), CD163 (Opal 620 Texas Red), and CD20 (Opal 488, FITC), with DAPI (blue) as nuclear counterstain. The tissue slides underwent a sequence of deparaffinization cycles (5 min each), then antigen retrieval was performed at 97 °C with using BOND Epitope Retrieval Solution 2 (Leica), primary antibody (45 min at 37 °C), incubation with multilinked HRPpolymer, incubation with fluorochrome-tyramid complex (Akoya Opal® Polaris 7-Color IHC Kit, NEL861001KT) (20 min at room temperature) and stripping of antibody complexes with ER2 solution (20 min at 97 °C), before the next immunolabeling cycle. After the last immuno-labeling cycle, spectral DAPI (Akoya Biosciences) was added to the slides for 5 min as a nuclear counterstain, and coverslips were mounted with ProLongTM Diamond Antifade Mountant (InvitrogenTM, ThermoFisher Scientific) mounting medium. All the steps previously described, except for the coverslips mounting, were performed on Leica Bond RX. Images were acquired using the Akoya PhenoImager HT 2.0 platform, and spectral unmixing was performed using the built-in synthetic spectra library to generate unmixed images for downstream analysis.

Digital image analysis was performed using QuPath (v0.5.1) [20] for single-cell segmentation and classification based on fluorescence intensity and localization. Based on the pathologist’s evaluation of H&E slide as a representative part of the whole tumor and multiple immunofluorescence analysis, an initial set of 25 candidate regions of interest (ROIs) per patient (500 μm in diameter) was manually selected to accurately target the invasive margin. Quality control filtering excluded regions with imaging artifacts or >60% DAPI-only cells. This threshold was applied to retain regions containing enough phenotyped cells for reliable spatial neighborhood analysis, while excluding areas that could not be robustly characterized with the current marker panel. In total, 198 ROIs were analyzed across the cohort corresponding to 15–21 representative ROIs per patient. Individual cells were segmented using StarDist algorithm [21] and annotated based on lineage-specific marker expression after supervised thresholding and classifier tuning. Across all ROIs, 397,270 single cells were segmented and annotated. Data were exported at single-cell and ROI levels for downstream analysis in a custom Python pipeline. Further details on the computational pipeline are provided in the Data and code availability statement.

For each patient, we pooled single-cell phenotype labels across all ROIs and computed the relative fraction of each cell phenotype (CD4⁺, CD8⁺, CD20⁺, CD163⁺ and DAPI-only). Immune cell fractions were retained as patient features for clustering analysis. To stabilize variance and normalize feature scales, fractions were logit-transformed and standardized to zero mean and unit variance. Patients were then grouped by unsupervised agglomerative clustering using Euclidean distance and Ward’s linkage (cophenetic correlation = 0.617). The optimal number of clusters was determined by dendrogram analysis. Given the non-normal distribution of cell fractions (Shapiro-Wilk test p < 0.05), differences between clusters were assessed using Mann-Whitney U tests with Benjamini-Hochberg correction and significance at p < 0.05.

Beyond composition, we explored the spatial organization of the TIME by characterizing local cellular microenvironments through cellular neighborhood (CN) analysis. Each cell was assigned to a neighborhood defined by the composition of its 10 nearest neighbors. Local environments were clustered with the k-means++ algorithm. After testing multiple values of k, k = 8 was selected as the biologically optimal configuration, yielding eight CNs (CN0–CN7) with distinct immune cell enrichments. For spatial visualization of CNs within ROIs, Voronoi diagrams were constructed using methods adapted from [22]. Due to non-normal CN distributions (Shapiro-Wilk p < 0.05), inter-cluster differences were tested using Mann-Whitney U tests and Benjamini-Hochberg correction (p < 0.05).

We investigated enrichment and co-localization patterns between immune lineages. First, we assessed the local co-occurrence of immune cell types through a neighborhood enrichment analysis [23], measuring, for each immune lineage serving as reference, the enrichment of each phenotype around its cells. Each lineage was considered in turn as a reference to capture pairwise interaction patterns across all cell types. Enrichment scores were computed and plotted for each cluster separately, using squidpy’s nhood_enrichment() functions, which evaluate statistical significance through permutation-based testing, setting a 50 μm radius and n = 1,000 iterations.

Finally, proximity analysis was conducted using a custom-developed Python algorithm to quantify spatial co-occurrence patterns between immune cell types within variable radii (from 10 to 200 μm at 10 μm intervals). For each ROI, the algorithm computed pairwise distance matrices between all cells using Euclidean coordinates. Cell pairs within each distance threshold were identified, and co-occurrence frequencies were calculated as the proportion of each cell type found in proximity to every other cell type, normalized by row to obtain relative fractions. ROI-level proximity data were then aggregated by computing mean and standard deviation across ROIs within each patient cluster.

We blindly reviewed the H&E slides of the tumor to describe the immune microenvironment at the invasive margin, subdividing between those that showed lymphoid aggregates in proximity of the tumor and those that presented only loose inflammatory infiltration. The inflammatory population in the portal tracts in the background liver parenchyma was not considered for the analysis.

Categorical variables were summarized descriptively, and comparisons between groups were performed using Fisher exact test or Chi-square test as specified in Supplementary Table S1. Survival outcomes were estimated using the Kaplan-Meier method, with OS and PFS compared between groups using the log-rank test. Median survival and 95% confidence intervals were reported, and hazard ratios (HRs) were estimated using univariate Cox proportional hazards models.

Given the limited sample size and number of events, all survival analyses were considered exploratory and hypothesis-generating.

To profile the spatial immune architecture of BTC prior to systemic therapy, we analyzed tumor samples from 11 treatment-naïve patients (4 iCCA, 5 pCCA, 1 dCCA, and 1 GBC). Baseline demographic and clinical characteristics, including sex, age, TNM classification, and mutational status, are reported in Table 1.

The analytical workflow, from image acquisition to spatial modeling, is illustrated in Figure 1A (see Methods for details). Briefly, sequential H&E-stained sections were reviewed by a pathologist to classify tumors by growth pattern, yielding four mass-forming (Patients 2, 6, 7, and 10) and seven periductal-infiltrating cases (Patients 1, 3, 4, 5, 8, 9, and 11), and to identify the invasive tumor margins (Supplementary Figure 1A).

Following this morphological guidance, whole-slide mIF scans were acquired and processed for fluorescence alignment. To capture the heterogeneity of the tumor-immune interface, we selected 25 representative ROIs per patient at the invasive margin. After quality filtering, 198 ROIs were retained, and approximately 400,000 individual cells were annotated. Representative tissue sections and ROIs are shown in Figure 1B; Supplementary Figure 2A.

Downstream analyses were performed to quantify immune lineage abundance, perform unsupervised clustering, and reconstruct spatial cellular neighborhoods at the invasive margin.

We first quantified the relative abundance of CD4⁺ T cells, CD8⁺ T cells, CD20⁺ B cells, CD163⁺ macrophages across patients, revealing marked inter-patient variability in immune composition (Figure 2A). To investigate whether patients could be grouped according to their immune profiles, we applied bottom-up hierarchical clustering to immune cell type fractions using Euclidean distance and Ward linkage (Figure 2B). This analysis separated the cohort into two major clusters: Cluster 1 (CL-1, n = 6 patients) and Cluster 2 (CL-2, n = 5 patients). We then compared immune composition between clusters. CL-1 was enriched in CD4⁺ T cells and CD163⁺ macrophages (Mann-Whitney U test, p < 0.0001 and p < 0.0001, respectively), whereas CL-2 was characterized by higher proportions of CD8⁺ T cells and CD20⁺ B cells (p = 0.025 and p < 0.0001, respectively) (Figure 2C).

We next asked whether compositional differences were paralleled by distinct spatial organizations. To test this, we applied cellular neighborhood (CN) analysis, in which each cell was assigned to a neighborhood defined by the phenotypic composition of its ten nearest neighbors. Across all samples, we identified eight discrete CNs (CN0–CN7), each representing distinct patterns of local cell interactions (Figure 2D, left). The defining characteristics of these niches included: B-cell dominance (CN1 and CN2, with CN1 also showing a CD4⁺/CD8⁺ T-cell component), CD4⁺ T-cell enrichment (CN4 and CN6), and macrophage or CD8⁺ T-cell prevalence (CN5 and CN7, respectively), while CN0 and CN3 displayed more heterogeneous compositions without a dominant lineage. A representative ROI is shown with Voronoi tessellation illustrating CN assignment (Figure 2D, right).

We then compared CN distribution across clusters. Patients in CL-2 exhibited a higher prevalence of CD20⁺ B cells that frequently co-localized with CD8⁺ T cells (CN1 and CN2), suggestive of tertiary lymphoid structure (TLS)-like arrangements (Mann-Whitney U test, p < 0.0001, <0.0001 and 0.032, respectively). In contrast, CL-1 was enriched in regions dominated by CD4⁺ T cells and CD163⁺ macrophages with limited cytotoxic infiltration (CN4, CN5 and CN6; all p < 0.0001) (Figure 2E). We displayed CN assignments in two representative patients using Voronoi diagrams for each cluster: Patient 9 (CL-1) and Patient 5 (CL-2) (Figure 2F). Comprehensive CN maps for all other patients are shown in Supplementary Figure 3.

Although the comparison was performed on a limited number of cases, all tumors lacking TLS-like aggregates and showing only mild or loosely organized inflammatory infiltrates on H&E evaluation belonged to CL-1 (see representative example in Figure 2G). In contrast, tumors displaying an inflammatory infiltrate organized in TLS-like structures were predominantly associated with the CL-2 phenotype (4 of 5 cases; see representative example in Figure 2H).

To further examine spatial interactions between immune cells, we performed a neighborhood enrichment analysis [23]. In CL-1, B cell neighborhoods were enriched for both B cells and CD4⁺ T cells, with only limited representation of CD8⁺ T cells. By contrast, in CL-2, B cell neighborhoods retained enrichment in CD4⁺ T cells but showed reduced self-enrichment of B cells (Figure 3A). Notably, co-occurrence of CD8⁺ T cells and CD163⁺ macrophages was observed in CL-1 but absent in CL-2 (Figure 3A).To investigate these relationships across broader spatial scales, we next performed proximity analysis. Relative fractions of neighboring cell types were quantified across increasing spatial radii and aggregated across ROIs within each cluster, considering all cell types included in the analysis. As shown in Figure 3B, immune spatial organization at the invasive margin diverged between clusters. In particular, in CL-1 the immune lineages CD8, CD20 and CD163 exhibited closer spatial proximity to CD4⁺ cells, suggesting a central role for CD4⁺ T cells in local immune organization. Conversely, in CL-2 cluster, CD20⁺ B cells emerged as the main spatial hub, with all other immune subsets showing preferential proximity to CD20⁺ cells. These patterns indicate distinct spatial organization across clusters, with differences in the relative proximity of CD4⁺, macrophage and B cell populations.

The study cohort included tumors arising from different anatomical sites (iCCA, pCCA, dCCA, and GBC), and cluster assignment was not associated with tumor location, with no statistically significant differences between groups (Supplementary Table 1).

We then evaluated potential associations between these immune spatial patterns and clinical or molecular tumor features. Targeted NGS was successfully performed on a subset of cases with available tumor material (n = 9/11 patients), revealing a limited number of alterations across a panel of cancer-related genes (Table 1). Within this subset, among CL-2 patients with available molecular data (3/5), all cases harbored co-occurring TP53 and ARID1A alterations (3/3), whereas no such co-occurrence was observed in CL-1 (0/6) (p = 0.0119, Fisher exact test) (Figure 2A).

We next evaluated whether cluster assignment was associated with clinical outcome (Table 2). Kaplan–Meier analysis did not reveal statistically significant differences in OS or PFS between clusters (Figures 4A,B, left), with log-rank p values of 0.82 and 0.76, respectively. Descriptively, median OS was 24.1 months in CL-1 and 13.5 months in CL-2, while median PFS was 3.2 months and 5.0 months, respectively. However, these differences were not supported by the overall shape of the survival curves. Consistently, univariate Cox models did not identify a significant association between cluster assignment and survival outcomes (OS HR = 1.201; PFS HR = 0.796).

Additional univariate Cox regression analyses for OS and PFS, were performed across clinical variables, namely growth pattern (periductal vs. mass-forming), grading (G2 vs. G3) and treatment regimen (CG vs. CGD) (Figures 4A,B, right); none showed statistically significant associations with OS or PFS.

To explore treatment-associated correlates of outcome in relation to immune phenotype, we stratified the cohort according to overall survival above or below the cohort median (23.9 months). Within the CG subgroup, this yielded numerically balanced long- and short-survival groups, allowing exploratory comparison of baseline immune features (Figures 5A,B). Long survivors showed higher fractions of CD163⁺ macrophages and lower CD20⁺ B cells at baseline, and were enriched in regions characterized by high CD163⁺ macrophage content with limited lymphoid representation (CN5) and depleted in B-cell–rich regions (CN2). The same analysis was not performed in the CGD subgroup because only one patient fell into the short-survival category.