
==== Front
iScience
iScience
iScience
2589-0042
Elsevier

S2589-0042(24)01893-5
10.1016/j.isci.2024.110668
110668
Article
Multiplex imaging of localized prostate tumors reveals altered spatial organization of AR-positive cells in the microenvironment
Ak Çiğdem 1211
Sayar Zeynep 12811
Thibault Guillaume 29
Burlingame Erik A. 210
Kuykendall M.J. 1
Eng Jennifer 2
Chitsazan Alex 1
Chin Koei 1
Adey Andrew C. 13
Boniface Christopher 1
Spellman Paul T. 13
Thomas George V. 14
Kopp Ryan P. 15
Demir Emek 17
Chang Young Hwan 2
Stavrinides Vasilis 6
Eksi Sebnem Ece eksi@ohsu.edu
1212∗
1 Cancer Early Detection Advanced Research (CEDAR), Knight Cancer Institute, OHSU, Portland, OR 97239, USA
2 Department of Biomedical Engineering, School of Medicine, OHSU, Portland, OR 97209, USA
3 Department of Molecular and Medical Genetics, Knight Cancer Institute, OHSU, Portland, OR 97239, USA
4 Department of Pathology & Laboratory Medicine, School of Medicine, OHSU, Portland, OR 97239, USA
5 Department of Urology, School of Medicine, Knight Cancer Institute, Portland, OR 97239, USA
6 Cancer Institute, University College London, London, UK
7 Division of Oncological Sciences, School of Medicine, OHSU, Portland, OR 97239, USA
∗ Corresponding author eksi@ohsu.edu
8 Present address: Zeynep Sayar, Parse Biosciences, Seattle, WA, 98109 USA

9 Present address: Guillaume Thibault, Merck, 126 East Lincoln Avenue PO Box 2000. Rahway, NJ 07065 USA

10 Present address: Erik Burlingame, Johnson & Johnson Innovative Medicine, New Brunswick, NJ, 08933 USA

11 These authors contributed equally

12 Lead contact

03 8 2024
20 9 2024
03 8 2024
27 9 1106681 4 2024
19 7 2024
1 8 2024
© 2024 The Authors
2024
https://creativecommons.org/licenses/by-nc-nd/4.0/ This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).
Summary

Mapping the spatial interactions of cancer, immune, and stromal cell states presents novel opportunities for patient stratification and for advancing immunotherapy. While single-cell studies revealed significant molecular heterogeneity in prostate cancer cells, the impact of spatial stromal cell heterogeneity remains poorly understood. Here, we used cyclic immunofluorescent imaging on whole-tissue sections to uncover novel spatial associations between cancer and stromal cells in low- and high-grade prostate tumors and tumor-adjacent normal tissues. Our results provide a spatial map of single cells and recurrent cellular neighborhoods in the prostate tumor microenvironment of treatment-naive patients. We report unique populations of mast cells that show distinct spatial associations with M2 macrophages and regulatory T cells. Our results show disease-specific neighborhoods that are primarily driven by androgen receptor-positive (AR+) stromal cells and identify inflammatory gene networks active in AR+ prostate stroma.

Graphical abstract

Highlights

• We provide a detailed spatial map of single cells via cyclic imaging

• Distinct mast cell populations are spatially associated with M2 macrophages and Tregs

• Grade-specific cellular neighborhoods are identified based on stromal AR expression

• Stromal AR expression is correlated with NF-κB1 and HIF1α

Microenvironment; Cancer; Omics

Subject areas

Microenvironment
Cancer
Omics
Published: August 3, 2024
==== Body
pmcIntroduction

Prostate tumors are composed of epithelial cells that exhibit molecular and cellular differences and interact with a complex ecosystem of stromal cells: fibroblasts, immune cells, mesenchymal stem cells, smooth muscle, blood vessels, and innervating neurons.1 The complexity and heterogeneity of the disease are indicated by pathological and single-cell omic evidence.2,3,4 It is essential to simultaneously dissect this cellular heterogeneity and generate a spatial map between cell types and states to better stratify patients into treatment groups and to identify more effective immunotherapy strategies. Not surprisingly, these cell states, cell types, and bidirectional cellular interactions may not be solely recovered from hematoxylin and eosin (H&E) staining, even with advanced machine learning algorithms. Single-cell omic approaches can identify cell types and states, but not the spatially coordinated behavior of heterogeneous cells. Immunohistochemistry is gaining widespread adoption in clinical practice to mitigate this problem, but it is often deployed for a few markers in a non-systematic, non-standardized manner.

Recent advancements in multiplex proteomic profiling and image analysis facilitated the profiling of >30 proteins from thousands of single cells on single tissue sections, opening a new path toward predicting disease subtypes based on immune microenvironments, linked to progression.5,6 In breast cancer, structured immune microenvironments linked to survival were identified.7,8,9 In colorectal cancer, distinct cellular communities, which can distinguish molecular subtypes and aid in the risk stratification of aggressive forms of cancer, were revealed.10 These studies show the power of mapping spatial associations in the tumor microenvironment (TME) for predicting patient trajectory and revealing novel targets for therapeutic interventions.4

In the prostate TME, dynamic changes in the physical microenvironment of prostate tumors, such as hypoxia, wound healing, and chronic inflammation, are postulated to contribute to cancer initiation and progression.11,12,13 These processes require coordination between immune and non-immune stromal cells. However, our understanding of non-immune stromal cell types is limited, especially around androgen receptor (AR) expression in the prostate TME. A few studies report loss of stromal AR expression that is strongly correlated with the worst patient outcomes.14,15 However, a better delineation of grade-specific AR-positive (AR+) stromal cell types and of pathways involved in AR regulation in stromal cells is needed.16

Here, we developed an extensive antibody library to probe changes in the prostate TME and identified 27 epithelial, immune, and non-immune stromal cell phenotypes that exhibit prostate cancer grade-specific differences in spatial composition (Figure 1A). We used cyclic immunofluorescent (cycIF) imaging to spatially map 699,461 cells from 13 patients and three clinical groups (Figure 1A). Further, we analyzed copy-number aberrations (CNAs) for 10 patients to investigate changes in the spatial coordination of immune and non-immune stromal cells. We show AR+ stroma-driven neighborhood organization associated with clinical groups and define transcriptional networks active in AR+ stroma using single-cell combinatorial indexing assay for transposase-accessible chromatin sequencing (sci-ATAC-seq) from the same group of patients.17 Our study provides a comprehensive spatial and regulatory map of the prostate TME for patients with low and high-risk localized prostate tumors and highlights the importance of spatial profiling in identifying patient-specific disease subtypes.Figure 1 Single-cell phenotypes of primary prostate tumors identified via cyclic immunofluorescent imaging

(A) Diagram showing the experimental (top) and computational (bottom) pipeline of cyclic immunofluorescent (cycIF) analysis of primary prostate tumors (image created with Biorender.com). Tissue sections were stained, imaged, and quenched repetitively for 10 rounds. Cell segmentation detected 699,461 cells in total that were annotated into 27 cell types following a multi-step clustering approach. 10 cellular neighborhoods, which were composed of distinct recurrent cell types, were identified across all clinical groups.

(B) The antibody library used in the study distinguishes between epithelial (red), immune (yellow), and non-immune stromal cells (teal) and includes several functional status markers (pink) for each cell type.

(C) Example cycIF images from TAN, G3, and G4 tissue sections show αSMA (blue), AR (red), and Ki67 (magenta) markers and H&E images from the adjacent tissue section. White arrows show normal (TAN) or cancerous (G3 and G4) glands. Yellow arrows show stroma.

(D) CellPose segmentation algorithm determines the boundaries of each nucleus (left) and cell (right). Arrows point out to the defined boundaries of a stromal cell.

(E) Numbers of cells were segmented from each clinical group TAN, G3, and G4.

Results

Developing a prostate-cancer-specific antibody library to probe immune and non-immune stromal cell type changes

To study systemic changes in the prostate TME, we developed an extensive antibody library to capture prostate cancer cells in relation to immune and non-immune stromal cell types. We tested over 150 antibodies and validated more than 100 antibodies following a three-step verification for each antibody (Figure S1, Data S1 and S2). Our antibody library is specifically geared toward capturing changes in prostate cancer biology.

We selected 28 targets from this panel for cycIF, which allowed for a thorough classification of specific cell types in the prostate tissue: epithelial, immune, and non-immune stromal cells. We also included molecular and epigenetic markers, functional status markers, and a proliferation marker to capture specific cell states associated with patient outcomes (Figures 1B and S2).

Spatial profiling of 699,461 cells isolated from tumors spanning three clinical groups

Formalin-fixed paraffin-embedded whole-tissue sections from patients with localized primary prostate tumors were obtained for cycIF analysis (Figures 1A and 1C). Whole-slide images varied between 2 and 26 mm2 in surface area. We did not generate tumor microarrays (TMAs) to avoid any sampling bias.18 Tissue sections were assigned to clinical groups: tumor-adjacent normal (TAN), primary Gleason grade 3 (G3), or primary Gleason grade 4 (G4), based on their histopathological evaluations. If a tissue section contained >95% normal glands or stroma and <5% G3 glands, it was assigned to the TAN group. If >75% of the section consisted of G3 tumors, it was assigned to the G3 group, and if >75% of the section consisted of G4 tumors, it was assigned to the G4 group (Figures S3 and S4). The Cancer of the Prostate Risk Assessment Post-Surgical (CAPRA-S) scores for all patients were calculated based on clinical patient data.19 We identified five patients with low-, six patients with intermediate-, and two patients with high-risk CAPRA-S scores (Figures S3 and S4).

We analyzed the co-expression of 28 proteins across 13 whole-tissue sections using cycIF. Single-cell segmentation via CellPose identified 699,461 cells, with similar numbers of cells from each clinical group (G3, G4, and TAN) (Figures 1D and 1E).20

The density of epithelial, immune, and non-immune stromal cell types indicates changes in the cellular organization of tissues associated with clinical group

Our previous work uncovered a loss of epigenomic heterogeneity in epithelial prostate cancer cells.17 To further investigate cell type heterogeneity, we performed a multi-step clustering analysis for an in-depth classification of epithelial, immune, and non-immune stromal cells in prostate tumors. First, we clustered all 699,461 single cells based on the expression of all 28 markers in our panel and detected 49 clusters (Figure 2A).21 We used uniform manifold approximation and projection (UMAP) to visualize single cells based on Louvain clustering results, colored by either patient IDs (left) or cell types (right) (Figure 2A).22 Single-cell clusters on the UMAP space showed patient-specific epithelial cell clusters (Figure 2A). In contrast, immune and non-immune stromal cells across all patients formed a single, mixed cluster separate from all epithelial cells (Figure 2A).Figure 2 Combined expression of four markers, AR, AMACR, H3K27ac, and H3K4me3, distinguishes prostate epithelial cells across clinical groups

(A) Louvain clustering of 699,481 cells shows clinical-group (left) and cell-type-specific (right) clusters.

(B) The fraction of epithelial cells (purple) increased, while immune (yellow) and non-immune stromal (green) cells decreased in G4 tumors.

(C) Annotations of epithelial cell clusters based on the clinical groups (left) and cell types (right).

(D) Shannon entropy results measure the distribution of cells expressing epithelial markers across clinical groups, showing an overall decrease in marker heterogeneity with increasing primary tumor grade.

(E) The combined expression of four markers (AR, AMACR, H3K27ac, and H3K4me3) distinguishes between G4 samples vs. TAN and G3 samples (chi-square test with p < 2.2e−16). The mosaic plot shows deviations of observed frequencies as compared to expected. The area of the box is proportional to the difference in observed and expected frequencies. Pink indicates smaller and blue indicates larger observed frequencies as compared to expected (p value<2.22 × 10−16).

(F) A G4 grade gland detected by an H&E shows expression of AR, AMACR, H3K27ac, and H3K4me3, but does not express CK8 (luminal) or CK5 (basal) markers (left). A G3 grade and a normal gland detected by an H&E show CK8 and CK5 expression.

(G) Plot showing the spatial distribution of H3K27ac and AMACR in an example TAN (top) and G4 (bottom) samples.

To delineate cellular heterogeneity, we plotted Z-scored protein expression values as a heatmap and used epithelial (CK5, CK8, ECAD, and chromogranin A), immune (CD3, CD4, CD11b, CD20, CD45, CD66b, CD68, CD163, FOXP3, HLADRB1, and Tryp) and non-immune stromal (αSMA and VIM) markers to create three exclusive categories (Figure S5). Further details about cell phenotyping are provided in the STAR methods section. This analysis identified 435,702 epithelial cells (62% of cells), 209,521 stromal cells (30% of cells), and 54,238 (8% of cells) immune cells (Figure 2B). We analyzed the percentage of cell types across three clinical groups, which showed an increase in epithelial cells by 21% and 34% and a decrease in immune and non-immune stromal cells by 15% and 27% in G3 and G4 tumors, respectively, compared to TAN samples (Figure 2B).

The combined expression of epithelial markers shows epigenetic differences between clinical groups

We sought to further define cellular heterogeneity for epithelial, immune, and non-immune cells. First, we investigated epithelial cell heterogeneity, by re-clustering only the epithelial cells using the protein markers NCAM, chromogranin A, CD44, VIM, CK5, NKX3.1, AMACR, CK8, ECAD, H3K27ac, H3K4me3, and AR, and we identified five epithelial cell types: basal, luminal, AMACR-high luminal, neuroendocrine, and VIM+ epithelial cells that are likely in the process of epithelial-to-mesenchymal transition (EMT) (Figure 2C).

AMACR is consistently overexpressed in prostate cancer epithelium and is suggested as a cancer-risk-associated marker.23,24 We observed that 86% of all luminal cells express average levels of AMACR, and we found a separate cluster of luminal cells with high AMACR expression in G3 and G4 samples that we labeled as “AMACR-high” (Figure 2C).

While we observed an increase in the percentage of proliferating (Ki67+) AMACR-high luminal cells in G4 samples, the density of Ki67+ AMACR-high luminal cells was similar in G3 and G4 tumors (Figure S5). Significantly, we found an increase in the density of Ki67+ VIM+/EMT cells in G4 tumors, but not in TAN and G3 (Figure S5).

To quantify the heterogeneity of epithelial cells, we measured Shannon entropy for each clinical group.21 We observed a non-significant decrease in heterogeneity from TAN to G3 and G3 to G4, suggesting epithelial cells in TAN samples show more differences in epithelial marker expression, as opposed to G4 samples (Figure 2D).

Cell types are static categories that may not reflect the heterogeneity of cell states in prostate cancer epithelia. To determine the specific markers associated with the heterogeneity of epithelial cell states, we investigated the combinatorial expression of 12 epithelial markers in 699,461 single cells (Figure S6). Strikingly, we found that the combined expression of AMACR, H3K27ac, H3K4me3, and AR in epithelial cells showed a significant difference across the three clinical groups: 21.5% in G4, 6.4% in G3, and 0.4% in TAN (Figure S6). This significance was not observed for CK8+ luminal cells: approximately similar percentages of epithelial cells, 20.3% in G4, 18% in G3, and 19.8% cells were positive for CK8, AMACR, H3K27ac, H3K4me3, and AR (Figure S6). We performed a chi-squared test that showed that the statistical difference was significant between the three clinical groups (Figure 2E). We confirmed the loss of CK8 in G4 tumors using H&E images of individual TAN, G3, and G4 glands (Figure 2F).

H3K27ac marks active regulatory regions associated with increased transcription.25,26 We observed a 10% increase in the number of cells with H3K27ac expression in G4 tumors (Figure S6). To determine whether the spatial distribution of this marker may contribute to disease progression, we calculated Moran’s I values for all epithelial markers (Figure S6) and visualized spatial expression patterns on the tissues (Figure 2G). Moran’s I uses continuous marker expression across a tissue to determine spatial auto-correlation.27 We observe broad, low-level expression of H3K27ac that overlaps with AMACR expression in G4 samples, as opposed to TAN (Figure 2G).

Identification of Tryp+ mast cell populations that differentially express CD44, CD90, and GZMB

We further stratified immune cells and identified 13 myeloid and lymphoid cell types with Louvain clustering based on the expression of 17 immune markers: CD66b, CD11b, CD20, NCAM, PD1, CD163, HLA-DRB1, CD68, CD90, Tryp, CD44, GZMB, CD45, CD3, CD4, FOXP3, and AR (Figure 3A). Tryp+ mast cells constituted the second largest population of immune cells (18%) following helper T cells (19%), forming two distinct clusters on the UMAP, which suggested heterogeneous mast cell populations in localized prostate tumors (Figures 3A and 3B).Figure 3 Mast cells in the prostate tumor microenvironment differentially express CD44, CD90, and GZMB

(A) UMAP shows immune cell subtype clusters identified via Louvain clustering. Mast cell clusters 1 and 2 are shown in gray.

(B) Fractions of immune cell subtypes show high numbers of mast cells across all clinical groups.

(C) Heatmap showing log-normalized signal intensity values of each immune marker for the two mast cell clusters identified in our dataset. A large fraction of mast cells in cluster 2 (orange) are GZMB+ (blue).

(D) Reclustering of mast cells, UMAPs colored by clinical group (left), and previously identified mast cell clusters 1 (green) and 2 (orange) (right).

(E) The association plot reveals a significant abundance of mast cluster 2 in G3 samples and mast cluster 1 in TAN samples and no difference in G4. Rectangles in each row are positioned relative to a baseline indicating independence. If the observed frequency of a cell is greater than the expected one, the box rises above the baseline; otherwise, it falls below. The area of the box is proportional to the difference in observed and expected frequencies (p value <2.22 × 10−6).

(F) The scatterplot shows a positive correlation between GZMB and CD90 expression in mast cells (R = 0.29).

(G) Bar plot showing the fraction of GZMB+ and GZMB− cells across mast clusters.

(H) Boxplot showing AR, CD44, CD90, PD1, and CD68 expression differs significantly in mast clusters 1 and 2 (values represent mean ± SEM; t test ∗∗∗∗p < 0.001).

(I) Percentage distributions of mast cells showing mono-, bi-, and tri-expression of CD44, CD90, GZMB, and AR based on mast cell cluster IDs (top) and clinical groups (bottom).

(J) An example cycIF image that shows heterogeneous mast cell populations, AR+GZMB-CD90+CD44, AR-GZMB+CD90+CD44+, AR-GZMB-CD90+CD44+, AR+GZMB-CD90+CD44+. AR is shown in red, AMACR in white, tryptase (mast cells) in cyan, GZMB in purple, CD90 in green, and CD44 in yellow. The scale bar indicates a length of 100 μm.

Even though the presence of large fractions of Tryp+ mast cells is repeatedly reported in localized prostate tumors, our understanding of their role in mediating tumor response is limited.28 We determined that the majority of mast cells express CD44, a receptor of hyaluronic acid that is abundant in the extracellular matrix (Figure 3C).29 We investigated this further by analyzing the differential expression of immune markers in two mast cell populations (Figure 3D). Our results show a significant abundance of mast cluster 2 in G3 patients and mast cluster 1 in TAN patients and no difference in G4 patients (Figure 3E). Furthermore, we identified a positive correlation between CD90 and GZMB expression in mast cluster 2 (Figure 3F).

Mast clusters 1 and 2 showed a striking difference based on the CD90 and GZMB status (Figures 3G and 3H). GZMB induces anti-tumor responses and is traditionally found in effector T cells.30 However, certain GZMB+ cells, such as myeloid-derived suppressor cells and mesenchymal stem cells (MSCs), can secrete GZMB to degrade receptors on effector T cells and induce immunosuppressive pro-tumor effects.30 Mast cluster 2 exhibited a higher fraction of GZMB+ mast cells (55%) as opposed to cluster 1 (21%) (Figure 3G). We identified differential expression of CD68, AR, PD1, CD44, and CD90 in two clusters of mast cells (Figure 3H). We investigated single, dual, and triple positivity for markers CD44, CD90, GZMB, and AR in Tryp+ mast cells and plotted these mast populations across clinical groups (Figures 3I and S7). We validated the presence of dual and triple-positive mast cell populations in the prostate TME using patient cycIF images (Figure 3J).

GZMB− and AR+ mast cells showed spatial association with M2 macrophages and Tregs in the prostate TME

To investigate whether identified mast cells are associated with pro- vs. anti-tumor immune cell types, we examined spatial co-occurrence patterns of mast cells in the prostate TME (Figure 4A). We divided mast cells into two groups based on their GZMB expression and analyzed the frequency of spatial co-occurrence of GZMB+ or GZMB− mast cells with all identified epithelial, immune, and non-immune stromal cell types.10 GZMB+ and GZMB− mast cells did not show any specific enrichment in their spatial auto-correlation in TAN samples (Figures 4B and S8). However, we identified an enriched spatial co-occurrence between GZMB+ mast cells and M2 macrophages and between GZMB− mast cells and CD90+ MSCs in G3 tumors (Figure S8). GZMB+ mast cells similarly showed highly significant spatial association with M2 macrophages and granulocytes in G4 tumors (Figures 4C and S8). We did not observe GZMB− cell spatial association with MSCs, but with neuronal cells in G4 prostate TME (Figure 4C).Figure 4 Mast cell subtypes spatially associate with distinct immune cells

(A) Voronoi diagrams of cellular neighborhoods for patients TAN_1, G3_3, and G4_1. Each dot shows mast cell distribution across whole tissue sections.

(B) Average Moran’s I score of immune cell subtypes across three clinical groups does not show differences based on mast cell marker expression (Tryp) (bottom, orange).

(C) Heatmaps of likelihood ratios of direct cell-cell contacts between 27 cell types in G4 tumors. Mast cells (Tryp+) were divided into two groups GZMB+/− (left) and AR+/− (right). GZMB+ mast cells showed highly enriched spatial correlation with M2 macrophages and are closely associated with granulocytes. AR+ mast cells showed highly enriched spatial association with Treg cells.

(D) Ripley’s L results of patient G3_3 highlight cell types that are spatially auto-correlated, such as CD4+CD3+ T cells (red), and others that are spatially dispersed, such as granulocytes (brown) and GZMB+ mast cells (black).

(E) Co-occurrence probability of each cell type given the presence of GZMB+ (left) and AR+ mast cells (right). The score is computed across increasing distances around each cell in the tissue. GZMB+ mast cells, non-APC M2 macrophage, AR+ mast cells, and CD4− CD3+ T cells show co-occurrence in close proximity.

(F) cycIF images showing GZMB+ mast cells (gray boxes) near M2 macrophages and GZMB− mast cells (white box). DAPI is shown in white, TRYP in blue, CD163 in yellow, and GZMB in red. The scale bar indicates a length of 50 μm.

Functional interactions of mast cells with regulatory T cells (Tregs) in colorectal cancer are shown to generate potently immunosuppressive and proinflammatory Tregs.31,32 To test whether expression profiles of mast cells impact their spatial correlation, we compared the behavior of mast cells with and without AR expression in the prostate TME. We identified AR+ mast cells that closely interact with Tregs in G3 and G4 prostate tumors as compared to TAN samples (Figures 4C and S8). Additionally, AR+ mast cells had enriched spatial association with M2 macrophages in G3, but not in G4 tumors (Figures 4C and S8). AR− mast cells did not show any specific enrichment in their spatial association (Figure S8). Taken together, our results show changes in the spatial association of mast cells with specific lymphoid, myeloid, and stromal cells based on their GZMB and AR expression, suggesting possible functional roles for distinct mast cell states.

To investigate whether mast cell spatial auto-correlation directly influences these results, we first generated Voronoi diagrams of each tissue section and then mapped mast cells onto Voronoi neighborhoods (Figure 4A).10 Next, we computed Moran’s I values for each patient, indicating the spatial auto-correlation of continuous variables, using mean signal intensity values of all markers in the panel (Figure 4B). Our results did not show any differences in mast cell auto-correlation associated with a clinical group (Moran’s I average = 0.08 for all groups). This provides further evidence that the specific spatial association of mast cells with M2 macrophages is not an artifact of spatial auto-correlation but is directly linked to tumor progression in high-grade tumors (Figures 4B and 4C).

We validated this result with Ripley’s L statistics, which is similar to Moran’s I but employs discrete values as its input, such as cell types.33 Mast cells did not show spatial auto-correlation, whereas CD4+ T cells did (Figure 4D). Additionally, we calculated the co-occurrence probability score for each cell type in comparison to GZMB+ (Figure 4E left) and AR+ (Figure 4E right) mast cell populations. GZMB+ mast cells exhibited significant spatial co-occurrence with non-APC M2 macrophages, and AR+ mast cells showed significant spatial co-occurrence with CD3+CD4− cytotoxic T cells based on increasing distance around each cell in the tissue (Figure 4E). We observed that mast cells close to CD163+ M2 macrophages frequently expressed GZMB (GZMB+Tryp+; white box on the left), while mast cells spatially distant from M2 macrophages did not express GZMB (GZMB-Tryp+; gray box on the right) (Figure 4F).

AR+ non-immune stromal cells drive cellular neighborhood organization in tumor samples

We further explored potential predictive factors linked to clinical groups by focusing on non-immune stromal cells. We delineated heterogeneous stromal cell populations by clustering them in isolation based on the expression of the following markers: αSMA (smooth muscle), VIM (mesenchymal fibroblasts and endothelial cells), CD90 and CD44 (MSCs and cancer-associated fibroblasts), NCAM (neurons), HLA-DRB1 (infiltrating antigen-presenting cells), and AR (Figures 5A and 5B).34,35Figure 5 AR-positive non-immune stromal cells direct cellular neighborhood organization in prostate tumors

(A) UMAP shows non-immune stromal cells colored based on clinical grade specific (left) and subtype (right).

(B) Fractions of non-immune stromal cells across clinical groups. TAN samples show high numbers of AR+ stromal cells as compared to G3 and G4 tumors.

(C) Recurrent cell neighborhoods (rCNs) are identified based on 27 cell subtypes. The top four neighborhoods are defined by an increased occurrence of AR+ stromal cells, the bottom six neighborhoods are defined by a lack of AR+ stroma.

(D) The frequency of each neighborhood is analyzed per clinical group. rCN1 and rCN2 (AR+ stromal cells, neuroendocrine and immune subtypes) frequency is significantly higher in TAN as compared to G4 and G3 tumors (t test p value 0.017 and 0.016), respectively. CN6 (AMACR-high luminal/other APC immune cells) frequency is significantly lower in TAN as compared to G4 tumors (t test ∗∗p value 0.001). See Figure S13 for cell type differential enrichment within rCNs. Scatterplot showing the frequency of each neighborhood per patient (bottom). Each colored dot represents a patient. rCN1, rCN2 ,and rCN6 show significant differences between TAN and G3 samples; rCN1 and rCN6 show significant differences between TAN and G4 (t test ∗p value<0.05).

(E) Dot plot of p values for the rCN frequency comparing prostate tumors across 13 iterations of individual patient exclusion (horizontal red line shows p = 0.05).

(F) Scatterplots of p values for the rCN frequency comparing prostate tumors across 100 iterations of random subsamples of 50% of the spatial data (horizontal red line shows p = 0.05).

(G) Bar plots showing the percentages of the Voronoi area of each neighborhood cover in TAN, G3, and G4 samples (right) and the actual Voronoi area of each ten CNs per grade (left).

(H) Voronoi diagrams of G4_1 and G4_5 patients showing rCN distributions. G4_1 and G4_5 patients show enrichment for very different rCNs, luminal (blue) vs. AR+ stroma (orange, green, brown, and yellow), respectively.

Our results uncover three populations of AR+ stromal cells: AR+ smooth muscle, AR+ αSMA, and VIM+ and AR+ other stromal cells (Figures 5A and 5B). We determined a significant increase in AR+ stromal cells (27%) in TAN samples, as compared to G3 and G4 tumors (6%) (Figure 5B). Interestingly, despite this increase, we did not observe significant enrichment in the spatial auto-correlation of AR+ stromal cell markers in TAN samples, whereas AR+ stromal cells in tumor samples showed significant enrichment in their spatial auto-correlation (Figure S9).

To determine the impact of AR+ stromal cell populations on tissue organization, we mapped the spatial frequencies of all 27 cell types for each patient and used non-hierarchical clustering to identify recurrent cell neighborhoods (rCNs)10 (Figure 5C). Each rCN was then mapped back onto the tissue sections to visualize cell type and neighborhood compositions.10

Our results indeed show that AR expression in prostate TME strongly informs the 10 rCNs identified across all patient samples. Four neighborhoods were dominated by AR+ non-immune stromal cells interacting with neuroendocrine (rCN1), mesenchymal/endothelial (rCN8), immune (rCN2), and smooth muscle and M2 macrophages (rCN5). The other six neighborhoods were dominated by the lack of AR+ stroma: luminal epithelial (rCN0), EMT/infiltrating VIM+ (rCN9), AMACR-high luminal cells interacting with APC immune cells (rCN6 and rCN7), basal and neuroendocrine (rCN3), and basal, neuroendocrine, and luminal cells (rCN4) (Figure 5C).

We sought to determine if the significant heterogeneity of cell types impacted the identified neighborhoods. We calculated the frequency of each rCN for TAN, G3, and G4 samples. Our results identified rCN1 and rCN2, which showed a statistically significant increase in TAN samples (Figure 5D). Both rCN1 and rCN2 are defined by recurrent spatial correlation of AR+ stroma with immune cells (Figure 5C). The frequency of these neighborhoods was significantly lower in both G3 and G4 samples as compared to TAN (Figure 5D). We observed that rCN5 was characterized by a notable lack of B cells, differentiating it from other stromal neighborhoods (Figure 5C). We also identified a significant increase in the frequency of rCN6, which consists of AMACR-high luminal cells in G3 and G4 tumors (Figures 5C and 5D). We repeated the neighborhood identification analysis 13 times, leaving one patient out each time (Figure 5E). Our results robustly captured the same rCN composition in each model: AR+ stroma-driven TAN vs. G3 and G4 differences were identified in all analyses, and AMACR-high-driven rCN change was identified in 11 out of 13 (Figure 5E).

To address the effects of spatial variability and sample size, we performed cellular neighborhood analysis on in silico tissues derived from each patient sample, with 10 or 50 tiles per clinical group (Figure S10). Our findings revealed consistent neighborhoods for TAN, G3, and G4 tissues, highlighting a significantly higher AR+ rCN frequency in TAN compared to G3 and G4 tumors (p value<0.05) (Figure S11). Additionally, we performed sub-sampling analysis through 100 iterations at four independent values (50%, 75%, 90%, and 99% of the data). The mean AR+ rCN frequency difference between groups (TAN vs. tumors) across these iterations was statistically significant (horizontal red line shows p = 0.05). In other words, even with 50% data subsampling, there is a statistically significant difference for AR+ stromal neighborhoods in TAN vs. G3 (86/100 of the iterations) and TAN vs. G4 (90/100 of the iterations) (Figure 5F). Combined, these results underscore the significance of AR regulation in non-immune stroma as a key factor in prostate cancer progression.

Loss of heterogeneity in cellular organization

Our results showed that AR+ stromal neighborhoods may distinguish between three clinical groups of prostate samples (Figures 5C–5F).12,36,37,38 To further delineate these clinical groups, we identified rCNs separately for TAN, G3, and G4 and identified cell types that are differentially enriched within specific neighborhoods across different clinical groups (Figures S8 and S10).

Looking at the frequency of each rCN in three clinical groups, we found TAN samples consist of neighborhoods that have similar frequencies across the tissue section (less than 25% frequency for each neighborhood) (Figure S9). In contrast, G4 samples consisted of one or two dominant neighborhoods (Figures 5D and S9). In G3 tumors, the predominant neighborhoods exhibited higher than 25%, but less than 40% frequency (Figures S9 and S12). Overall, these results suggest a loss of rCN heterogeneity associated with higher primary Gleason grade, specifically as a result of the decrease in AR+ stromal cell neighborhoods (Figures 5D and S8).

We calculated the percentage size distribution of each rCN per clinical group and showed a decrease in the actual size of the AR+ stromal neighborhoods in G3 and G4 tumors (Figure 5G). The size of the AMACR-high neighborhoods showed an increase in G3 and G4 patients as compared to TANs, supporting our results (Figure 5G). Distinct rCNs dominated patients with high-grade prostate tumors, suggesting specific molecular changes underlying the rCN organization (Figure 5H).

Identification of patient-specific genomic events

To further delineate patient-specific disease subtypes, we performed low-pass whole-genome sequencing for 10 patients and identified common CNAs. G4_1 and G4_2 patients did not have any known CNAs. G4_3 patient exhibited three known amplification events: NCOA2, NBN, and EGFR. G4_4 patient showed three known deletion events: PTEN, CDH1, and NKX3.1 (Figures S3 and S4). PTEN loss is associated with an immunosuppressive microenvironment that is mediated by Tregs and M2 macrophages,39 which correlates with our neighborhood analysis that shows a smaller number of APC M2 macrophages in G4 tumor neighborhoods as compared to TAN and G3 (Figure S13). Finally, G4_5 patient showed multiple amplification (PIK3CA and EGFR) and deletion (PTEN, CHD1, RB1, FOXO1, TP53, CDKN1B, and NKX3.1) events.

Overall, our combined genomic profiling and cycIF imaging results further stratify G4 tumor subtypes and exhibit major differences in TME organization that may be associated with known CNAs in prostate cancer, underlining the need for understanding patient-specific trajectories for precise treatment approaches (Figure 5H).

AR+ stromal cells in localized prostate tumors show activated NF-κB signaling

To gain further understanding of the biological targets and pathways in AR+ stromal cells, we shifted from imaging and focused on identifying potential gene networks active in AR+ cells in prostate TME. We used sci-ATAC-seq data from flash-frozen tumors collected from the same patients17 and identified myeloid, lymphoid, and non-immune stromal cells based on cisTopic clustering and gene expression scores40,41 (Figure 6A). We identified differentially accessible regions of AR+ vs. AR-null (AR0) cells and observed significantly enriched accessibility in chromatin regions that regulate nuclear factor κB (NF-κB) signaling in AR+ stromal cells (Figure 6B).42 Gene ontology (GO) analysis was performed using GREAT and the Panther classification system.42 To gain more mechanistic insights, we investigated transcription factor (TF)-binding motifs significantly enriched in AR+ vs. AR0 stromal cells using both Cistrome and JASPAR43 (Figure 6C). We identified significantly enriched TFs that are known AR regulators in AR+ cells such as SP1, YY1, and NFYA (Figure 6D). Interestingly, we observed that AR+ stromal cells in localized prostate tumors showed significant enrichment for the hypoxia marker HIF1α, suggesting a role for stromal AR regulation in hypoxia. We did not observe any enrichment for CXCR4, a known inflammatory marker, nor for FOXA1, a known repressor of HIF1α (Figure 6E). While there was no significant difference in IFNγ accessibility, a subgroup of AR0 cells was accessible for IFNγ, as expected (Figure 6D). Predicted gene scores confirmed that NF-κB1 was significantly more accessible in AR+ cells in the prostate TME (Figure 6D), validating our GREAT GO results (Figure 6B) and previous findings.2 We performed cycIF using AR, HIF1α, and NF-κB on three additional Gleason score 3 + 3 (G3 clinical group) patient samples and identified AR+ cells that also express HIF1α and NF-κB (Figure 6F). Overall, our sci-ATAC-seq analysis identified active transcriptional networks in both AR+ and AR0 stromal cells and defined utilizing markers that can be used to further investigate AR+ stromal cells in cycIF imaging.Figure 6 AR+ stromal cells in localized prostate tumors show activation of NF-κB signaling and hypoxia

(A) AR+ myeloid, lymphoid, and non-immune stromal cells were identified in single-cell ATAC-seq data from the same patients. Differentially accessible regions (DARs) are identified in AR+ vs. AR0 cells.

(B) GREAT GO results show an enrichment for NIK/NF-κB signaling pathway regulators in AR+ stromal cells.

(C) TF motif-binding sites identified in AR+ stromal cells show an enrichment for known AR regulators.

(D) Gene scores of identified TFs and pathways in single-cell ATAC-seq data in AR+ vs. AR0 cells.

(E) Heatmap showing fold enrichment values for DARs in AR+ vs. AR0 stromal cells identifies significant HIF1α and NF-κB accessibility in AR+ cells.

(F) IF staining of a G3 patient sample showing co-expression of HIF1α (red), AR (magenta), and NF-κB (green) staining in the prostate stroma. DAPI-marked nuclei are shown in blue on the left. AR+ stromal cells that express HIF1α (A), HIF1α, NF-κB (B), and NF-κB (C) are marked. Cell “A” provides an example of nuclear colocalization of HIF1α and AR, whereas in cell “B,” HIF1α and NF-κB are colocalized in the cytoplasm in a cell marked by AR+ expression in the nucleus. Cell “C” shows expression of NF-κB in the cytoplasm in an AR+ cell. The scale bar indicates a length of 50 μm.

Discussion

In this work, we integrated multiple single-cell approaches with bulk genomic sequencing and pathology to investigate changes in the prostate TME using patient samples across three clinical groups. We present novel spatial patterns within localized prostate tumors that imply patient-specific spatial profiles even within tumors of seemingly identical pathological grades. While 13 specimens are not representative of all known molecular subtypes of localized prostate cancer, our results show striking differences in disease subtypes within a primary Gleason grade, which would not be detected based on an H&E image. Future studies will involve further delineation of disease subtypes by generating an increased number of multi-modal spatial maps, matching H&E and cycIF images, and examining their predictive and prognostic clinical utility.

Our multi-modal study provides a key resource for the prostate cancer community to better understand the spatial associations between cancer and stromal cells in primary G3 and G4 tumors. While direct data integration was not performed, our approach encompassed multiple experiments on patient samples, including single-cell assay for transposase-accessible chromatin using sequencing (ATAC-seq), cycIF, H&E, and genomic DNA sequencing, improving our understanding of patient-specific differences prior to treatment.

Our multi-modal analyses underline the need for generating more detailed disease subtypes for patients with localized prostate cancer to improve prognosis. Biological underpinnings measured for each Gleason score 4 + 4 tumor exhibit differences in CNAs, immune response, and neighborhood organization. For instance, while we observed an overall increase in the frequency of immune-cold luminal neighborhoods in primary G4 tumors, we identified a patient dominated by an M2 macrophage-AR+ stromal cell neighborhood, which was consistent with previous studies showing T cell exclusion in tumors due to defects in antigen presentation, allowing tumors to evade the immune system.44,45 We identified one patient dominated by an M2 macrophage-AR+ stromal cell neighborhood. Interestingly, this patient was the only patient in our cohort with an RB1 deletion, which is a critical event associated with neuroendocrine prostate cancer.46 Although implications of AR+ stroma and immunosuppressive TME interactions need to be further characterized in the prostate, it has been reported that AR expression is correlated with M2 macrophage polarization and disease severity in infections such as asthma in the lung.47

Immunotherapies have been largely ineffective for prostate cancer.48,49 Prior studies have evaluated the enrichment of infiltrating mast cells, T cells, natural killer cells, and macrophages as predictors of prostate cancer outcome and better stratification of patients into accurate treatment groups.50,51,52 Although studies have consistently reported the prostate immune cell microenvironment as a largely immunosuppressive one,53,54,55 there are conflicting reports on the enrichment of specific cell types, such as mast cells.28,56,57,58 We believe that the inconsistencies between these reports are due to the lack of spatial information and evaluation of each cell type independently as opposed to the neighborhoods they are found within.54,59,60,61,62,63,64 A recent study,53 however, has taken spatial information into account by considering the compartment in which the cells are found (epithelia vs. stroma). Major strengths of this study are the use of multiplex digital pathology to quantify seven immune cell types and their large cohort size consisting of TMAs from over 900 patients. Comparison of immune cell density within the tumor epithelia vs. stroma showed a specific enrichment of Tregs and M1 and M2 macrophages in epithelial tumor regions. Consistent with other studies on Tregs45,55 and M2 macrophages,65,66 Andersen et al. also found that the enrichment of stromal Tregs and epithelial M2 macrophages were independent predictors of biochemical recurrence (BCR).53 Mast cells alone were not significantly associated with BCR.53 In our study, we used whole-tissue sections instead of TMAs to capture long-distance associations across 27 cell types and avoid sampling bias associated with the TMAs.18,67 We observed changes in the spatial association of mast cell populations with M2 macrophages and Tregs based on their GZMB and AR status. Our study revealed specific associations between AR+ mast cells with Tregs in G3 and G4 tumors and with M2 macrophages in G3 tumors. Furthermore, we identified a specific association of M2 macrophages with GZMB+ mast cells in G3 and G4 tumors. Overall, taken into consideration with the aforementioned studies, interactions of AR or GZMB-expressing mast cells with the immunosuppressive M2 macrophages and Tregs may be functionally important in promoting a more aggressive, BCR-prone prostate cancer phenotype. Functional interactions between mast cells and Tregs and macrophages mediated by cytokines have been investigated in mouse models of colorectal cancer.31,32 Other studies have shown interactions between M2 macrophages and CD90+ cells to promote a pro-tumor, invasive phenotype in prostate tumors.68,69 The interactions between mast cells and Tregs lead to potently immune-suppressive and proinflammatory Tregs, possibly through a positive feedback loop mediated through cytokines.70,71 Mouse models of prostate cancer cells do not have abundant mast cell populations, but we show these mast cell populations may play a role in immunosuppressive behavior in human adenocarcinoma. Co-culture studies are needed to investigate the functional interactions between GZMB+ and AR+ mast cell populations with M2 macrophages.

Mast cells are a suggested independent prognostic variable for prostate cancer progression.58 Several studies report associations between low tumor mast cell counts with poor outcomes such as BCR and metastasis.28,56,57,58,72 Although not statistically significant, our data also show a lower mast cell density in TAN compared to G3 or G4 samples. CD44 expression is shown to be a critical factor for mast cell proliferation, differentiation, and maturation in other diseases.73,74 Here, we report large fractions of mast cells expressing CD44, with no grade-specific changes in localized prostate tumors. In the context of myeloid malignancies, CD44 has been demonstrated as a mast cell invasion receptor and correlates with disease aggressiveness and tumor formation.74 It is possible that CD44 expression in mast cells modulates mast cell attachment to the extracellular matrix and contributes to extracellular matrix reorganization in high-grade tumors.29 Detailed markers are needed to identify CD44 isoforms that may be heterogeneously active across clinical groups. Overall, our analysis emphasizes the importance of evaluating the spatial context and molecular expression of cells as opposed to only cell counts when analyzing the TME, and we propose that the lack of consensus on which cell types are associated with adverse outcomes is due to the lack of spatial context in previous studies.

The impact of androgen deprivation therapy (ADT) is primarily considered through its effects on cancer cell proliferation and survival. Recently, the influence of ADT on T cell infiltration and function has gained attention.75 However, our understanding of AR regulation in other immune and non-immune stromal cells in prostate tumors, and how ADT may alter stromal response to cancer through AR regulation, remains limited. Mast cells in the skin have been shown to express AR in large fractions, with studies indicating that testosterone treatment does not impact degranulation.76 In cutaneous neurofibroma, activated AR promoted mast cell infiltration, and AR inhibition reduced mast cell infiltration.77 In prostate cancer, enzalutamide similarly increased mast cell infiltration by regulating AR expression, which was linked to neuroendocrine differentiation.78 These findings highlight the importance of investigating the impact of ADT on AR+ mast cells and their interactions with Tregs and M2 macrophages in creating a protumor microenvironment.

We also identified similarities between patients for each clinical group that may be used for automated patient stratification in future studies. We show that the combined expression of AMACR, H3K27ac, H3K4me3, and AR in epithelial cells can distinguish between TAN, G3, and G4, independent of CK8 status. Furthermore, we report a loss of heterogeneity in high-grade prostate tumors (G4 clinical group) as compared to low-grade (G3) and TAN samples. Our first line of evidence for loss of heterogeneity comes from our previous single-cell epigenomic analyses that demonstrate increased similarity in chromatin accessibility profiles of high-grade tumors.17,26,79 Here, we observed increased epithelial expression of epigenetic markers H3K27ac and H3K4me3, which are associated with activated transcription and accessible chromatin, accompanied by a decrease in Shannon entropy in high-grade tumors. Additionally, using cycIF and spatial statistics, we observed the specific areas of increased similarity in cancer cells on patient sections, identifying overlapping expression of H3K27ac and AMACR in G4 tumors.67,80 Finally, we report patient-specific rCNs in high-grade prostate tumors, showing that decreased heterogeneity is also reflected at the tissue-organization level.

Loss of stromal AR in advanced prostate cancer is associated with biochemical relapse and poor prognosis.15,16,36,81 However, there are no previous reports of delineating non-immune stromal cell types in human adenocarcinoma based on AR expression. While we detected 10 distinct stromal cell populations in this study, we observed many unclassified non-immune stromal cells that were marked positive for AR. We report that this “AR+ other stroma” cell population played a major role in neighborhood organization across three clinical groups. We suspect there will be tremendous value in stratifying these “AR+ other stroma” cells to better understand their role in neighborhood organization in low- vs. high-risk patients and in response to AR inhibitors. In future studies, additional markers will be added to cycIF panels to specifically stratify AR+ stromal cells, such as markers identified through sci-ATAC-seq.

sci-ATAC-seq provided valuable insights into active pathways in AR+ stromal cells. The number of cells captured via sequencing remains low as compared to the number of single cells captured via imaging. However, the depth of cis- and trans-regulator information gained from sci-ATAC-seq is unmatched in multiplex imaging, which is limited to the number of protein markers that can be probed from a single tissue section. Integrating single-cell sequencing and imaging approaches allowed us to overcome the limitations of each experimental tool and unlocked mechanistic insights related to NF-κB and hypoxia activation in AR+ stromal cells. Our findings validate the previous bulk RNA sequencing,82 proteomics,83 and single-cell RNA sequencing2 results that show an increased NF-κB signaling in immune cells isolated from localized prostate tumors that are AR+. We validated the co-expression of AR and NF-κB in prostate tumor stroma from patient samples via immunofluorescent imaging, which further supports the robustness of our findings. Furthermore, this suggests that, despite variations in tissue preparation methods between single-cell ATAC-seq and cycIF imaging, our results reveal consistent gene signatures that are preserved across different modalities. Recent advancements in high-throughput single-cell omics technologies will allow researchers to understand the role of stromal AR expression in regulating inflammation.

Limitations of the study

While our study presents compelling evidence regarding the role of AR+ stroma in prostate tumor organization supported by robust statistical analyses, several limitations should be acknowledged.

Our patient cohort consists of 13 individuals, which may appear limited. Although this sample size aligns with similar studies in spatial imaging and single-cell analysis,2,3,84,85,86 larger cohorts could provide more generalizable results. We performed extensive power analyses, including in silico tissue generation, leave-one-patient-out analysis, subsampling, and permutation tests to mitigate concerns about sample size and variability. These analyses demonstrated the robustness and stability of our findings. However, future studies with larger cohorts could potentially reveal more specific associations between AR+ stromal changes and genomic aberrations.

The leave-one-patient-out analysis was used to assess the impact of individual variability, which resulted in 11 out of 13 replications, in which the role of AR+ stroma in prostate tumor organization was highly significant. However, there is still a possibility that specific, unmeasured patient characteristics could influence these results. Our study attempted to mitigate this through rigorous statistical methods, but inherent biological variability remains a challenge.

Spatial heterogeneity within tissue samples poses another significant challenge. Accurate assessment of disease-relevant tumor structures requires the statistical power of whole-slide imaging (WSI), as opposed to TMAs, which consist of a small subsection of specimens.18,67 The requirement for WSI in a research and diagnostic setting comes with substantial cost. Per-patient datasets are over 100-fold larger than those with TMAs, and cohorts are more difficult to acquire since whole blocks must be accessed and re-cut. In contrast, 2D WSI (100 K cells per specimen) largely overcomes this problem (number of effective cells > 100) for the characterization of local neighborhoods. WSI captures the most variability to demonstrate the differences in distinct grading and their cell compositions. WSI is also the standard in conventional pathology and is regarded by the Food and Drug Administration as a diagnostic necessity.67

In summary, while our study provides significant insights into the role of AR+ stroma in prostate tumor organization and employs rigorous statistical analyses to support these findings, limitations related to sample size, generalizability, and methodological approaches should be considered when interpreting the results. Future studies with larger and more diverse cohorts, along with continued methodological advancements, could further validate and expand upon our findings.

STAR★Methods

Key resources table

REAGENT or RESOURCE	SOURCE	IDENTIFIER	
Antibodies	
	
alphaSMA (1A4)	Santa Cruz Biotechnology	Cat# sc-32251; RRID:AB_262054	
AMACR (Polyclonal)	Abcam	Cat# ab212378; RRID:AB_3096027	
AR_(CST) (D6F11)	Cell Signaling Technology	Cat# 8956S; RRID:AB_11129223	
CD11b_555 (EPR1344)	Abcam	Cat# ab206616; RRID:AB_2728809	
CD163_647 (EPR14643-36)	Abcam	Cat# ab218294; RRID:AB_2943126	
CD20 (EP459Y)	Abcam	Cat# ab198943; RRID:AB_2905499	
CD3_750 (EP449E)	Abcam	Cat# ab213608; RRID:AB_3096030	
CD4_647 (EPR6855)	Abcam	Cat# ab196147; RRID:AB_2923526	
CD44_555 (EPR1013Y)	Abcam	Cat# ab216647; RRID:AB_3096031	
CD45_750 (EP322Y)	Abcam	Cat# ab214437; RRID:AB_3096032	
CD66b_750 (G10F5)	Fisher Scientific	Cat# NBP2-80664; RRID:AB_3096017	
CD68_750 (KP1)	BioLegend	Cat# 916104; RRID:AB_2616797	
CD90 (EPR3133)	Abcam	Cat# ab216449; RRID:AB_2889264	
ChromA_750 (EP1030Y)	Abcam	Cat# ab215276; RRID:AB_3096034	
CK5_488 (EP1601Y)	Abcam	Cat# ab193894; RRID:AB_2893023	
CK8 (EP1628Y)	Abcam	Cat# ab192468; RRID:AB_2890258	
ECAD (EP700Y)	Abcam	Cat# ab201499; RRID:AB_2801591	
FOXP3 (PCH101)	Thermo Scientific	Cat# 14-4776-82; RRID:AB_467554	
GZMB 750 (EPR20129-217)	Abcam	Cat# ab219803; RRID:AB_2910576	
H3K27ac (EP16602)	Abcam	Cat# ab245911; RRID:AB_3096035	
H3K4me3 (C42D8)	Cell Signaling Technology	Cat# 11960S; RRID:AB_2797779	
HIF-1 alpha (ESEE122)	Novus Biologicals	Cat#NB100-131; RRID:AB_350068	
HLA DRB1 (EPR6148)	Abcam	Cat# ab207066; RRID:AB_3096026	
Ki67_CST (D3B5)	Cell Signaling Technology	Cat# 12075S; RRID:AB_3096029	
NCAM (EP2567Y)	Abcam	Cat# ab215981; RRID:AB_3096036	
NF-κB (D14E12)	Cell Signaling Technology	Cat#8242S	
NKX3.1 (3 9)	Thermo Fisher Scientific	Cat# 35–9700; RRID:AB_2533232	
PD1 (EPR4877(2))	Abcam	Cat# ab201825; RRID:AB_2728811	
Tryp (AA1)	Abcam	Cat# ab2378; RRID:AB_303023	
VIM (EPR3776)	Abcam	Cat# ab194719; RRID:AB_2909595	
	
Biological samples	
	
Tissue sections from patients	OHSU’s Biolibrary	IRB#4918	
	
Chemicals, peptides, and recombinant proteins	
	
AlexaFluor (AF) 488, AF555, AF647 and AF750 dyes	ThermoFisher Scientific	A20006	
AlexaFluor NHS ester reagents	ThermoFisher Scientific	A20006 and A20000	
Amicon Ultra spin columns	Sigma-Aldrich	UFC5010	
Ampure XP Beads	Beckman Coulter	A63880	
Citrate, pH 6	Sigma-Aldrich	C-1909	
IGEPAL CA-630	Sigma-Aldrich	I8896-100ML	
Target Retrieval Solution, pH 9	Agilent	S2367	
PBS, pH 7.4	Fisher Scientific	BP39920	
Normal Goat Serum (NGS)	Fisher Scientific	S-1000	
Bovine Serum Albumin (BSA)	ThermoFisher Scientific	AM2616	
Slowfade Gold DAPI	ThermoFisher Scientific	S36938	
#1.5 thickness coverslip	Corning Life Sciences	2980-243 and 2980-24	
Protease inhibitor tablet	Roche	Cat. 11873580001	
Kapa HiFi HotStart ReadyMix	Fisher Scientific	KK2602	
Bioanalyzer HS dsDNA Kit	Agilent	5067–4626	
QiaQuick PCR purification kit	Qiagen	28104	
	
Deposited data	
	
Raw and analyzed data	This paper	Figshare+ DOI: https://figshare.com/s/116766067fa229670e17	
	
Software and algorithms	
	
Image Preprocessing – conversion to TIFF	Zeiss	Zen Blue Microscopy Software	
Image Preprocessing	Quantitative Imaging Systems, LLC	QiTissue software	
Nuclei and Cell Segmentation	CellPose	https://github.com/MouseLand/cellpose	
Signal Intensity Normalization	RESTORE	https://gitlab.com/Chang_Lab/cycif_int_norm	
Cell typing with GPU accelerated PhenoGraph	Grapheno	https://gitlab.com/eburling/grapheno	
Heatmap visualization		https://gitlab.com/eburling/BCTMA/-/blob/master/pyviz_figs.py	
Cluster Validation		https://github.com/zeynepsayar/PCa_immune/blob/main/mapping.m	
Tissue Area Calculation	NIH, Bethesda, Maryland	ImageJ software	
Cellular Neighborhood	Schürch et al.10	https://github.com/nolanlab/NeighborhoodCoordination	
squidpy	Palla et al.33	https://github.com/scverse/squidpy	
BWA-MEM		https://github.com/lh3/bwa	
IchorCNA		https://github.com/broadinstitute/ichorCNA	
cistrome		http://cistrome.org/	
chromVAR		https://bioconductor.org/packages/release/bioc/html/chromVAR.html	
motifmatchr		https://bioconductor.org/packages/release/bioc/html/motifmatchr.html	
Creating In Silico Tissues and Power Analysis	This paper	Github: https://github.com/cigdemak/cycIF_prostate
Zenodo: https://doi.org/10.5281/zenodo.10333970	
snapATAC	Fang et al.41	https://github.com/r3fang/SnapATAC	
	
Other	
	
Single-cell ATAC sequencing data	Eksi et al.17	GEO: GSE171559	
JASPAR database		https://jaspar2020.genereg.net/	
Cistrome	Liu et al.81	http://cistrome.org/	
TintoRetriever Pressure Cooker	BioSB		
Axio Scan.Z1 slide scanner	Zeiss		
Colibri7 LED light source	Zeiss		
Orca Flash 4.0 V3 CMOS camera	Hamamatsu Photonics		
Plan-apochromat 20x/0.8 objective	Zeiss		
Sonicator for DNA fragmentation Covaris E220	Covaris Inc		
2100 Bioanalzyser system	Agilent Technologies		
NextSeq 500 Sequencer	Illumina		
BD FacsAria Fusion cell sorter FACSDiva v8.0.3	BD		
QuantStudio v1.7.1 for RT-PCR	ThermoFisher Scientific		

Resource availability

Lead contact

Further information and requests for resources and reagents should be directed to and will be fulfilled by the lead contact, Ece Eksi (eksi@ohsu.edu).

Materials availability

This study did not generate any new reagents.

Data and code availability

• Imaging data has been deposited to Figshare+ (https://doi.org/10.25452/figshare.plus.25199468) and is publicly available as the date of the publication. All cycIF data deposited follows the multiplex imaging data standards published by the Human Tumor Atlas Network (HTAN): https://humantumoratlas.org/standard/imaging. This paper analyzes existing, publicly available single-cell ATAC sequencing data (GSE171559). The accession number for the dataset is listed in the key resources table. Any additional information required to analyze the data reported in this work is available from the lead contact upon request.

• Custom codes have been deposited on Github and Zenodo (Github: https://github.com/cigdemak/cycIF_prostate, https://github.com/zeynepsayar/PCa_immune; Zenodo: https://doi.org/10.5281/zenodo.10333970) and are publicly available as of the date of the publication.

• Any additional information required to reanalyze the data reported in this work is available from the lead contact upon request.

Experimental model and study participant details

Human subjects

Thirteen white male patients, aged between 46 and 81, with treatment-naïve, localized prostate cancer were identified from the OHSU’s Biolibrary using medical records and reviews of H&E images. Tissue samples were obtained from men undergoing radical prostatectomy with informed consent through Biolibrary under IRB#4918. Participants provided samples voluntarily without compensation. Our research adheres to all pertinent ethical guidelines set forth by the Knight Cancer Institute at OHSU. Five unstained 5-micron FFPE adjacent tissue sections and an adjacent H&E stained tissue section from each patient were obtained. H&E images were reviewed by a pathologist and a Gleason score was assigned to each patient. Patients were classified into risk groups based on The Cancer of the Prostate Risk Assessment Post-Surgical (CAPRA-S) scoring system (Figures S3 and S4).87

Method details

Cyclic IF antibody panel

Antibodies used in the study can be found under the Antibodies section in the key resources table and are also listed in Data S1. Validation of antibodies for cycIF to mark cell types in the prostate tumor microenvironment is shown step by step in Figure S1. Tissue sections used in antibody validations are shown in Data S2. Antibodies were either obtained as commercially conjugated to fluorophores or in a BSA and Azide free buffer and conjugated to AlexaFluor (AF) 488, AF555, AF647, and AF750 dyes in-house (ThermoFisher Scientific, A20006) [21].

Fluorophore conjugation to antibodies

Antibodies were either obtained as commercially conjugated to fluorophores or in a BSA and Azide free buffer and conjugated to AlexaFluor (AF) 488, AF555, AF647, and AF750 dyes in-house. AlexaFluor NHS ester reagents (ThermoFisher Scientific, A20006, and A20000) were dissolved in DMSO to a final concentration of 10 mM. Antibodies were reacted with fluorophore dyes at a 1:10 antibody-to-dye molar ratio in 0.1 M sodium bicarbonate buffer of pH 8.3. Conjugation reactions were carried out on a rocker for 2 h at room temperature in the dark. A buffer exchange step using 10 kDa Amicon Ultra spin columns (Sigma-Aldrich, UFC5010) was performed to remove the excess Alexa Fluor dye from conjugation reactions.

FFPE tissue processing

FFPE tissue sections were deparaffinized as follows17: The paraffin embedding was removed from the tissue sections via an overnight incubation at 50°C followed by a 1-h incubation at 65°C. Slides were immediately transferred into Xylene and sequentially immersed in a fresh Xylene solution two times, 5 min each; 100% ethanol two times 5 min each; 95% ethanol two times 2 min each; 70% ethanol two times 2 min each and distilled water two times 5 min each. Antigen retrieval was done in a medical decloaking chamber (BioSB, TintoRetriever Pressure Cooker). Tissue slides were placed in a Coplin jar containing 10 mM citrate buffer, pH 6 (Sigma-Aldrich, C-1909), and incubated at high pressure for 15 min. Each slide was then immersed in hot distilled water and transferred into 1× Target Retrieval Solution, pH 9 (Agilent, S2367) for 15 min. Following antigen retrieval, tissues were briefly washed in ddH2O, followed by 5 min in 1× phosphate-buffered saline (PBS), pH 7.4 (Fisher, BP39920).

Cyclic IF imaging of prostate tissues

Cyclic IF protocol was adapted with minor revisions88,89 as follows: Before immunofluorescence staining, the slides were incubated in a quenching solution (10% 10× PBS, 0.4% 5M NaOH, 3% H2O2) under broad spectrum light for 1 h to reduce inherent tissue autofluorescence. After quenching, tissues were blocked with 5% NGS (Fisher Scientific, S-1000) in 1% BSA (Thermo Scientific, AM2616) solution for 30 min at room temperature. For each round of cyclic IF staining, primary antibodies were diluted in 5% NGS in 1% BSA. Slides were stained overnight at 4°C in a humid chamber. Antibodies were washed four times in 1× PBS for 5 min. Tissue sections were coverslipped using Slowfade Gold DAPI mounting media using a #1.5 thickness coverslip (Corning Life Sciences, 2980-243 and 2980-245) and imaged using a Zeiss Axio Scan.Z1 slide scanner. After a successful scan was obtained, the fluorophore signal was quenched as described above under broad-spectrum light for 1 h. Slides were washed 3 times in 1× PBS for 5 min. Quenching was confirmed under a microscope. Primary antibody staining, imaging, and quenching steps were repeated for 10 rounds, using up to four distinct antibodies conjugated to AF488, AF555, AF647, and AF750 and DAPI nuclear stain in each round. The order in which each antibody was used in the cyclic IF experiment can be found in Figure S2B. Quenched signal images were taken following rounds 3 and 10 and before round 1.

Whole slide scanning

An Axio Scan.Z1 microscope (Zeiss, Germany) with a Colibri7 LED light source with Cy7 line (Zeiss, Germany) was used for automated high-throughput whole slide scanning. LED intensities were kept constant for all rounds of imaging, where the LED exposure times were adjusted in each round. Exposure times were kept constant between all tissues in each round. Filter sets for DAPI, GFP, Cy3, Cy5, and Cy7 were used for 5-channel imaging in each round of cyclic IF. Images were taken with an Orca Flash 4.0 V3 CMOS camera (Hamamatsu Photonics, Japan) at 20× resolution using a plan-apochromat 20×/0.8 objective (Zeiss, Germany). Image acquisition profiles were set using the ZEN Blue edition imaging software, where the coarse and fine focus strategies and light exposure times were manually adjusted for each round of cyclic IF imaging.

Low-pass WGS library preparation and sequencing

Low-pass WGS libraries were created using the Kapa Biosystems Hyper Prep kit (KAPA Biosystems, Capetown, South Africa) using 11.8 ng/mL of G4_1; 17.6 ng/mL of G4_3; 13.8 ng/mL of G4_4; G4_4 10.8 ng/mL; 13 ng/mL G4_5; 13.3 ng/mL G3_1; 9.2 ng/mL G3_2; 11.4 ng/mL G3_3; 9.6; G3_5; 11 ng/mL TAN_1; 10.1 ng/mL TAN_2; 10.4 ng/mL TAN_3 DNA as input. DNA was fragmented by sonication to 150 bp using a Covaris E220 (Covaris Inc., Woburn, MA, USA) and ligated to dual-indexed Illumina-compatible sequencing adaptors using KAPA Hyper Prep Reagents and protocol and then PCR amplified for 9 cycles using Illumina library amplification primers. Amplified libraries were assessed for size distribution and concentration using the 2100 Bioanalyzer system (Agilent Technologies, Santa Clara, CA, USA) with HS dsDNA Kits. Samples were then sequenced the Illumina NextSeq 500, paired-end 75 bp with dual 14-bp indexing cycles. Raw fastq files were demultiplexed using in-house software, aligned using BWA-MEM, and copy number alterations were called using the IchorCNA software package using default parameters90 (https://github.com/broadinstitute/ichorCNA).

Quantification and statistical analysis

Single-cell ATAC sequencing analysis

We used the single-cell ATAC sequencing data from 13 matching patients over 18 patients that was published previously.17 Radical prostatectomy samples from 13 patients with low- and high-risk prostate cancer were obtained as flash-frozen samples and were prepared for single-cell ATAC sequencing using a combinatorial indexing platform.17,91 Sequencing base calls were converted to fastq files using bcl2fastq. Reads were pre-processed using snaptools snap-pre command with default parameters and aligned to the hg19 genome using snaptools align paired-end command using default parameters and bwa specified as the aligner.41 Peaks were called using runMACS command from snapATAC R package with the following MACS2 parameters: -nomodel -shift 100 -ext 200 -qval 5e-2 -B –SPMR.41 Clusters consisting of immune or non-immune stromal cells were identified based on cluster IDs and gene expression scores.17,41 We included all subjects and data without exclusion, maintained blinding to all patient information, and determined sample sizes based on the availability of patient samples.

Motif analysis

We used the function “getJasparMotifs” from R package chromVar which fetches motifs from the JASPAR database.92 The function “matchMotifs” from the motifmatchr R package finds which peaks contain which motifs. We found enriched motifs of both AR+ and AR0 cell groups using their DARs. One option is the p.cutoff for determining the stringency of motif calling. The default value is set to 0.00005, which gives reasonable numbers of motif matches for human motifs. The function can also return the number of motif matches per peak and the maximum motif score per peak. We calculated enrichment p-values comparing motif scores per peak between AR+ and AR0 cell groups using the Kolmogorov-Smirnov test implementation of the stats library in R.

Cyclic IF image analysis workflow

Image preprocessing

All images were acquired as.czi files from the Axio Scan Z1 microscope (Zeiss, Germany). ZEN Blue image processing software (Zeiss, Germany) was used to convert.czi images to 64-bit BigTiff images without any compression of the original grayscale data for each channel. The following steps; registration, segmentation, and autofluorescence subtraction, were performed.

A custom image registration algorithm was used to overlay DAPI images from each round of cyclic IF imaging. Registration quality was confirmed by visualizing on QiTissue software (Quantitative Imaging Systems, LLC, Pittsburgh, Pennsylvania) or the Napari application in Python.

Background signal removal

Images acquired after the chemical quenching step following the third round of cyclF experiment were selected as the autofluorescence image. The autofluorescence image was subtracted from each round. Following autofluorescence subtractions, a background smoothing operation was performed to smooth out staining gradients and uneven illumination patterns on the images. This was done by applying white-top-hat operation with the structuring element radius equal to the biggest nucleus radius within each tissue.

Image segmentation

Cellpose was used for the segmentation of individual nuclei.20 For nuclei segmentation, the z-projection image of all DAPI images from each round of cyclic IF imaging was used. Cytoplasm segmentation mask was created by dilating the nuclear segmentation mask. Average signal intensity values from individual nucleus, cell membrane, cytoplasm (nucleus minus cell membrane), and whole cell (nucleus plus cell membrane) compartments were calculated.

To eliminate parts of tissue that degraded and fell off the slides, we created a binary mask from the DAPI image from the last round of imaging and determined and used only the unique cell labels present in the last round DAPI mask.

Signal intensity normalization and thresholding

Signal intensities of each marker were normalized by background fluorescence intensity. A multi-step approach was followed for background fluorescence determination. First, we used an automated algorithm called RESTORE, which uses mutually exclusive markers, to calculate background fluorescence intensities.93 Cutoff value outputs from RESTORE were visually inspected for each marker in each tissue and manually adjusted. Cutoff values were then used to normalize signal intensities by dividing the raw fluorescence intensity values by the cutoff value. Therefore, values above “1” were “marker positive” as their signal intensity values were above the background signal intensity. Based on normalized intensity values, each cell in the dataset was labeled as proliferating (Ki67+) or non-proliferating (Ki67-) using a threshold of normalized Ki67 intensity >1. Similarly, each immune cell was identified as cytotoxic (GZMB+) or non-cytotoxic (GZMB-) and AR expressing (AR+) or not (AR-).

Clustering and cell type identification

Cell type identification analysis was performed using publicly available Python and R packages.21 Single-cell phenotypes were identified based on the normalized signal intensity values from our 28-plex cyclic IF data. To cluster single cells, we used an unsupervised Louvain community detection algorithm called PhenoGraph (k = 40, to detect rare cell phenotypes).94 Before the application of PhenoGraph to our data, we performed 99th percentile normalization and arcsin transformation (cofactor = 5). The GPU-accelerated implementation of PhenoGraph was implemented from https://gitlab.com/eburling/grapheno.21 For heatmap visualization of single-cell phenotypes and respective marker expression, “clustermap” function from the seaborn package in Python was used from previously developed algorithms available at https://gitlab.com/eburling/BCTMA/-/blob/master/pyviz_figs.py.

Following PhenoGraph, clusters were visually confirmed by mapping the location of cells from each cluster onto segmentation masks. We selected random samples to validate cluster annotations by doing side-by-side comparisons of cells mapped onto segmentation masks and IF images. MATLAB code for mapping cells onto segmentation masks is available on GitHub (https://github.com/zeynepsayar/PCa_immune/blob/main/mapping.m).

Cell type annotation

We selected 28 markers which would allow for the thorough classification of specific cell types in the prostate tissue: epithelial (CK5, CK8, ECAD, EPCAM, Chromogranin A (ChromA) and CD44), stromal (aSMA, VIM, CD90, NCAM) and immune cells (CD3, CD4, CD11b, CD20, CD45, CD66b, CD68, CD163, FOXP3, HLADRB1, NCAM, Tryp). We also included molecular and epigenetic markers (AR, p53, ERG, NKX 3.1, AMACR, PTEN, H3K27ac, H3K4me27), immune cell functional status markers (GZMB, PD1) and a proliferation marker (Ki67). We used Tryptase as our mast cell marker. Markers used for each cell type annotation are shown in Figure S2.

We adopted a 3-step clustering and annotation strategy to obtain the most granular clustering results and be able to detect rare cellular subtypes. In the first step, we performed clustering in Phenograph using 40 K-nearest-neighbors and normalized signal intensity values from all 28 markers listed above, identifying 49 clusters. The 49 clusters were annotated into 2 groups: epithelial, or non-epithelial based on cytokeratin (CK5, CK8) and cadherin (ECAD) expressions.

As the second step, we used Phenograph to cluster non-epithelial cells and annotate each cluster as either immune or non-immune stromal cells. For this, we used 50 K-nearest-neighbors and all markers in our panel except for epithelial markers (CK5, CK8, ChromA, ECAD, AMACR). We detected 30 clusters and annotated each as either immune or non-immune stroma. Stromal clusters were identified based on the expression of aSMA and VIM, and immune clusters based on the expression of various CD markers.

As the third step, we ran Phenograph separately on epithelial, immune, and non-immune stromal clusters and re-clustered them using 50 K-nearest-neighbors to identify specific cell subtypes within each major cell type. After running each clustering analysis, we grouped clusters with similar protein expression profiles under the same annotation identifying a total of 27 distinct cell types in the prostate microenvironment.

For re-clustering of epithelial cells, the following epithelial markers were used: H3K27ac, NKX3.1, NCAM, AMACR, CK8, ECAD, H3K4me3, AR, VIM, CK5, CD44 and ChromA (Figure 2C). For the re-clustering of non-immune stromal cells, the following epithelial markers were used: AR, aSMA, HLA-DRB1, CD90, VIM, CD44, and NCAM. We identified 23 clusters and classified them into 10 subtypes (Figure 5A). For re-clustering of immune cells, the following epithelial markers were used: CD66b, CD11b, CD20, PD1, CD163, PD1, HLA-DRB1, AR, CD68, CD90, Tryp, GZMB, CD44, FOXP3, CD4, CD45 and CD3. We identified 27 clusters and classified them into 12 subtypes (Figure 3C).

Clusters were visually confirmed by mapping the location of cells from each cluster onto segmentation masks. We selected random samples to validate cluster annotations by doing side-by-side comparisons of cells mapped onto segmentation masks and IF images. MATLAB code for mapping cells onto segmentation masks is available on GitHub (https://github.com/zeynepsayar/PCa_immune/blob/main/mapping.m).

UMAP visualization

To enable the 2-dimensional visualization of our high-plex feature-by-cell data, we used the GPU accelerated, RAPIDS library adaptation of UMAP on 99th percentile normalized and arcsin transformed (cofactor = 5) data. The default set of parameters was used for both visualization methods.

Quantification of epithelial marker heterogeneity using shannon entropy

Cells from each tissue were labeled as positive for CK8, CK5, AR, AMACR, H3K27ac, H3K4me3 and NKX3.1 if their normalized signal intensities were above 1. A Shannon entropy figure (Figures 2D and S6) was generated using UpSetPlot as implemented in the following algorithm: https://gitlab.com/eburling/BCTMA/-/blob/master/rapids015_figs.py. The heterogeneity of epithelial markers was computed by measuring the Shannon entropy (Figures 2D and S6A).21 Statistical differences between marker heterogeneity in TAN, G3, and G4 clinical groups were quantified using the pairwise_tukey function from pingouin.95

Combined marker expression analysis

To test the significance of the association between clinical groups and combined expression of epithelial markers (AR, AMACR, H3K27ac, H3K4me3) we used mosaic plots. Mosaic plots visualize the deviance of observed combined marker expression in each clinical group from the expected combined expression frequencies of an independence model. The mosaic plot in Figure 2E was generated with “mosaic” function from R’s “vcd” package using default parameters. The size of boxes in the mosaic plot is proportional to the difference between expected and observed frequencies. Boxes are colored in red if observed frequencies are smaller than expected, and blue if observed frequencies are greater than expected. We performed a chi-square test with the null hypothesis that grades and combined marker expressions are independent. P-values <0.001 reject the null hypothesis and indicate grades and combined marker expressions are dependent.

Mast cell differential marker expression analysis

Normalized mean intensity values of CD68, AR, PD1, CD44, and CD90 were plotted for each cell in Mast Cluster 1 and Mast Cluster 2 from all samples in the cohort. t-test was used to identify any statistically significant differences between individual and combined marker expressions between Mast cluster 1 vs. 2. p-values <0.05 were significant.

Cell density calculations

Cell counts were calculated for each cluster in epithelial, stromal, and immune populations. Tissue areas were calculated in mm2 using ImageJ software (NIH, Bethesda, Maryland). The whole tissue region was selected as the region of interest (ROI) and the area was calculated in mm2 using the measure option and the manually set scale. The fraction of Ki67-positive proliferating cells within each immune, epithelial, and stromal cell types was calculated and plotted in Figure 5E. Frequencies of each cell phenotype found in each tissue section were plotted in Figure 3B. The statistical significance of differences between cell densities was tested using a t-test (p < 0.05).

Spatial analysis

Pairwise cell-cell interactions

Pairwise cell-cell interactions were calculated using the R Shiny app “app_CRC_contacts.R” (https://github.com/nolanlab/NeighborhoodCoordination).96 Cell types were omitted if the number of unique adjacent cells for that cell type was less than 100 cells. Calculated values were mapped using the ComplexHeatmap package in R. To calculate pairwise cell-cell contacts, we calculated the likelihood ratios of the 27 cell types that we identified as follows:Loglikelihoodratio=log2(Observedfrequencyofco-ocurranceTheoreticalfrequencyofco-ocurrance)

Voronoi diagrams were generated through the following process in https://github.com/nolanlab/NeighborhoodCoordination: FCS files obtained from VorteX were exported and then processed using a custom Java algorithm. This algorithm was designed specifically to create Voronoi diagrams and derive contact matrices between individual cells.

Cell neighborhood identification

Cell neighborhoods were determined by using existing analysis platforms found in https://github.com/nolanlab/NeighborhoodCoordination. The ten nearest neighbors of each individual cell were determined. Each window was represented as the frequency of each cell type found within the ten nearest neighbors. In total we identified 27 cell types, therefore each window was a vector of length 27. To determine the cell neighborhoods, these windows were clustered using K-means clustering with k = 5 and k = 10 with Python’s scikit-learn library and “MiniBatchK-Means" function. We identified significantly enriched neighborhoods in each clinical group by performing t-tests (p < 0.05).

To validate the reproducibility of cell neighborhoods identified from the dataset and their enrichment in specific clinical groups, we iteratively performed cell neighborhood analysis by leaving out one patient at a time. For example, we identified cell neighborhoods from all cells in the dataset excluding cells from patient 1 and performed t-test to identify which cell neighborhoods were significantly enriched in TAN vs. G3 vs. G4. We repeated this analysis for each patient in the cohort.

Enrichment of cell types within cell neighborhoods

To calculate differential cell type enrichment in rCNs, differential enrichment analyses using linear models was conducted. The models were estimated using the following equation: Yn,c = β0 + β1X + β3Yc + e. In this equation, Yc represents the logarithm of the overall frequency of cell type c, X is an indicator variable for the patient group, Yn,c represents the logarithm of the frequency of cell type c in rCN n, βi are coefficients, and e represents Gaussian noise with a mean of zero. To avoid computational issues, a pseudocount of 1e-3 was added before taking logarithms. The estimation of these models was performed using the statsmodels Python package. The resulting coefficient estimates and p-values for β1 were extracted and visually presented (significant if p < 0.05).

Spatially variable markers with Moran’s I

Moran’s I evaluate spatial autocorrelation; whether a marker is spatially clustered, dispersed, or random based on the spatial expression of that marker:I=nW∑i=1n∑j=1nwi,jzizj∑i=1nzi2

where zi is the deviation of the feature from the mean (xi−1n∑i=1nxi), wi,j is the spatial weight between observations, n is the number of spatial units and W is the sum of all wi,j.

The Moran’s I global spatial auto-correlation statistics were calculated using Python’s Squidpy package implementation.33 We computed the Moran’s I score with “squidpy.gr.spatial_autocorr” and mode = 'moran'. We first generated a spatial graph with “squidpy.gr.spatial_neighbors” which provides a score on the degree of spatial variability of marker expression. The statistic as well as the p-value were computed for each marker, and FDR correction was performed (significant if p < 0.05). We visualized the top markers after applying Moran’s I spatial autocorrelation statistics.

Spatially variable cell types with Ripley’s L

Ripley’s spatial statistics is a family of spatial analysis methods used to describe whether points with discrete annotation (i.e., cell-types) in an area of interest follow random, dispersed or clustered patterns. To calculate Ripley’s L, we used “squidpy.gr.ripley” function with mode = ‘L’.

Cell-type co-occurrence ratio

The co-occurrence ratio of cell types provides a score on the co-occurrence of the cell type of interest across spatial dimensions. It is defined as follows :p(exp|cond)p(exp)

Where p(exp|cond) is the conditional probability of observing a cell-typeexp conditioned on the presence of a cell-type cond, whereas p(exp) is the probability of observing exp in the radius size of interest. The score was computed across increasing radii size around each cell in the tissue. We used Python’s “squidpy” library, “squidpy.gr.co_occurrence” function to calculate the cell-type co-occurrence ratio.33

Spatial heterogeneity and spatial power analysis

We performed several rigorous statistical analyses to ensure the sufficient spatial power and consistency of our results: (1) creating in silico tissue to augment sample size and to address sample variability, (2) leave-one-patient-out analysis to assess the impact of patient-to-patient variability on the results, (3) subsampling analysis to estimate the variability in our results and showed the stability of our findings, (4) permutation test to show further benchmarking and robustness of our findings. In all of the analyses, we considered a p-value of 0.05 or less to be statistically significant. All these different tests are different ways to address statistical power to determine the significance of the results.

Creating in situ tissues – ISTs to augment the data

Several essential parameters are outlined for planning a spatial omics study with adequate statistical power.67 The authors presented a flexible approach for generating in silico tissues (IST) and combined it with spatial profiling datasets to establish an exploratory computational framework for analyzing spatial power. Using their experimental design, to address effects of spatial variability and sample size, we divided each tumor sample into multiple regions, then performed cellular neighborhood analysis on these multiple regions, and showed the consistency of our results.

We divided each whole-tissue image into multiple regions of size 2K-by-2K pixels after testing different size parameters: 500, 1K, 2K, 3K, 4K, and 5K. First, we selected the most populated (number of cells in the grid) top 10 ISTs for each grade from 2K-by-2K grids ending up with 30 ISTs in total. The number of cells these 30 ISTs contain ranges from 1973 to 3457 cells. We obtained 22556, 28611, and 33440 cells for TAN, G3, and G4 grades respectively, totaling 84607 cells. As a second experiment, we also selected the most populated (number of cells in the grid) top 50 tiles for each grade from 2K-by-2K grids ending up with 150 IST in total. The number of cells these 150 ISTs contain ranges from 1460 to 3929 cells. We obtained 91117, 121734, and 124326 cells for TAN, G3, and G4 grades respectively, totaling 337177 cells.

Leave one patient out analysis

To further confirm the robustness of our findings, we performed an exclusion analysis on the spatial neighborhood data. Leave-one-patient-out analysis helps in understanding how sensitive the results are to the presence or absence of specific patients. It is particularly useful when dealing with a small dataset and can provide insights into the robustness and stability of the findings. It also allows us to assess the impact of patient-to-patient variability on the results.

For each round of analysis, we excluded one patient’s data from the imaging dataset. Then, performed rCN analysis, calculated the rCN frequencies for each patient, and compared AR + rCN frequency between TAN and tumors. We repeated this process iteratively for each patient, leaving one out at a time. Hence, we performed the analysis 13 times, each time excluding a different patient’s data. We evaluated the consistency of our results across the iterations. We observed that rCN patterns were robust and not overly influenced by specific patients. We also confirmed that the AR + rCN frequency between TAN and tumors (G3 and G4) remained statistically significant 10 out of 13 iterations for TAN vs. G3, and 11 out of 13 iterations for TAN vs. G4.

Sub-sampling analysis

We performed subsampling to validate that our sample size was sufficient to detect clinically relevant and statistically robust differences in our AR + rCN conclusion. First, 100 iterations of subsampling with different seeds were performed at four independent values; namely 50, 75, 90, and 99% of the data (sub-sampled from each grade data with given sub-sampling percentage). Then we run the rCN algorithm on each sub-sampled data to identify the AR + rCNs. We next tested whether the mean AR + rCN frequency was significantly different between TAN and tumors (G3 and G4) within each iteration of the sampled data. The distribution of p-values for the mean AR + rCN frequency differences between patient groups for each of these 100 iterations shows that even with 50% data subsampling, most iterations yield a p-value <0.05 (Figures 5F, S10, and S11). The proportion of statistically significant p-values is 86, 82, 83 and 91% (for TAN vs. G3); and 90, 93, 96, and 92% (for TAN vs. G4) for subsamples of 50, 75, 90, and 99%, respectively).

Spatial significance test 1

Spatial t-test with permutations: We used scipy.stats.ttest_ind() function in Python with exact permutations; permutations = np.inf during rCN frequency comparison. This permutes the patient’s grade assignment while comparing rCN frequencies with t-tests. We also separately tested equal variance across grades verified using Bartlett’s test. AR + rCN (CN1) frequency was significantly different between TAN and tumors (G3 and G4) with p-values of 0.018 and 0.036, respectively.

Spatial significance test 2

Spatial permutation test with sample mean difference statistic: We performed a permutation test using scipy.stats.permutation_test(). We defined the difference between the sample means as a test statistic to be used in scipy.stats.permutation_test(). We performed an exact test, n_resamples = np.inf. The permutation distribution of the test statistic is provided for further investigation, not only for AR + rCN but for all rCN frequency comparisons. AR + rCN (rCN1) frequency was significantly different between TAN and tumors (G3 and G4) with p-values of 0.036 and 0.036, respectively.

Spatial significance test 3

Spatial negative control by shuffling the rCN IDs: In this test, first, we calculated the ground truth neighborhood distribution matrix, as we described in cell neighborhood identification section above. Then, the rCN assignment labels on the tissue are shuffled to relabel tissue by rCNs, preserving the tissue structure. At each shuffle, the observed neighborhood frequency per grade was recalculated. We repeated the shuffling 100 times, and the proportion of statistically significant p-values is 8 and 29%. This result shows that the fact the AR + rCNs are differentiative of TAN vs. tumors was not an artifact.

Spatial significance test 4

Spatial negative control by both sub-sampling 99% of the data and shuffling the rCN IDs: We first subsampled our single-cell spatial data at 99%, 100 times. Then we identified the rCNs. Then repeated 3.a., shuffled the rCN IDs. At each shuffle of the 99% sub-sampled data, the observed neighborhood frequency per grade was recalculated. The proportion of statistically significant p values is 6 and 14%.

Other statistical tests

All statistical tests to infer p-value for significant differences (p < 0.05) in mean were performed using the t-test.

Supplemental information

Document S1. Figures S1–S13

Data S1. Comprehensive list of antibodies utilized in this study, related to STAR Methods–Cyclic IF Antibody Panel

Data S2. Catalog of tissue sections employed in antibody validations, related to STAR Methods–Cyclic IF Antibody Panel

Acknowledgments

We are grateful to the Biolibrary and the Histopathology Shared Resource (HSR) at OHSU for providing patient samples, specifically Aletha Letsch and Cheyenne Martin. The Biolibrary and HSR were supported by 10.13039/100000002 NIH grants P30 CA069533 and P30 CA069533 13S5 through the Knight Cancer Institute. Our research complies with all relevant ethical regulations of the Knight Cancer Institute at OHSU. Samples from patients undergoing radical prostatectomy were collected with informed consent through Biolibrary under IRB#4918. Participants did not receive compensation. We thank the Spellman lab and the Human Tumor Atlas Network (HTAN) for engaging in discussions. We thank Stefanie Kaech Petrie and her team at the Advanced Light Microscopy Core at OHSU for their technical help. We thank Lindsey Cauthen for her feedback on the manuscript and her revisions. We also thank Shelley Barton for providing feedback on the manuscript. This project was supported by funding (CEDAR3410918 ) from the 10.13039/100020527 Cancer Early Detection Advanced Research Center at Oregon Health & Science University, Knight Cancer Institute (S.E.E.).

Author contributions

Z.S. conducted cycIF experiments. C.A. and Z.S. analyzed the cycIF imaging data. C.A. and A.C. analyzed the single-cell ATAC-seq data. C.B. and S.E.E. performed CNA experiments and analysis. G.T. reviewed pathology slides. R.P.K. and V.S. reviewed the clinical patient data. Z.S., C.A., and S.E.E. designed the experiments and drew conclusions from the data. S.E.E. supervised all experiments and analyses. J.E., G.T., E.A.B., E.D., and Y.H.C. provided analysis pipelines for cycIF registration, segmentation, and clustering steps. A.A. provided experimental and computational pipelines for single-cell ATAC-seq. Z.S., C.A., and S.E.E. wrote the manuscript with input from all authors. S.E.E. reviewed and edited the manuscript.

Declaration of interests

The authors declare no competing interests.

Supplemental information can be found online at https://doi.org/10.1016/j.isci.2024.110668.
==== Refs
References

1 Lawson D.A. Kessenbrock K. Davis R.T. Pervolarakis N. Werb Z. Tumour heterogeneity and metastasis at single-cell resolution Nat. Cell Biol. 20 2018 1349 1360 10.1038/s41556-018-0236-7 30482943
2 Tuong Z.K. Loudon K.W. Berry B. Richoz N. Jones J. Tan X. Nguyen Q. George A. Hori S. Field S. Resolving the immune landscape of human prostate at a single-cell level in health and cancer Cell Rep. 37 2021 110132 10.1016/j.celrep.2021 34936871
3 Hirz T. Mei S. Sarkar H. Kfoury Y. Wu S. Verhoeven B.M. Subtelny A.O. Zlatev D.V. Wszolek M.W. Salari K. Dissecting the immune suppressive human prostate tumor microenvironment via integrated single-cell and spatial transcriptomic analyses Nat. Commun. 14 2023 663 10.1038/s41467-023-36325-2 36750562
4 Rozenblatt-Rosen O. Regev A. Oberdoerffer P. Nawy T. Hupalowska A. Rood J.E. Ashenberg O. Cerami E. Coffey R.J. Demir E. The Human Tumor Atlas Network: Charting Tumor Transitions across Space and Time at Single-Cell Resolution Cell 181 2020 236 249 10.1016/j.cell.2020.03.053 32302568
5 Keren L. Bosse M. Marquez D. Angoshtari R. Jain S. Varma S. Yang S.R. Kurian A. Van Valen D. West R. A Structured Tumor-Immune Microenvironment in Triple Negative Breast Cancer Revealed by Multiplexed Ion Beam Imaging Cell 174 2018 1373 1387 10.1016/j.cell.2018.08.039 30193111
6 Gustafson K.T. Sayar Z. Le H. Gustafson S.L. Gower A. Modestino A. Ibsen S. Heller M.J. Esener S. Eksi S.E. cyc-DEP: Cyclic immunofluorescence profiling of particles collected using dielectrophoresis Electrophoresis 43 2022 1784 1798 10.1002/elps.202200001 35753078
7 Bodenmiller B. Imaging mass cytometry and multiplatform genomics define the phenogenomic landscape of breast cancer Nature Cancer 1 2020 163 175 35122013
8 Danenberg E. Bardwell H. Zanotelli V.R.T. Provenzano E. Chin S.F. Rueda O.M. Green A. Rakha E. Aparicio S. Ellis I.O. Breast tumor microenvironment structures are associated with genomic features and clinical outcome Nat. Genet. 54 2022 660 669 10.1038/s41588-022-01041-y 35437329
9 Johnson B.E. Creason A.L. Stommel J.M. Keck J.M. Parmar S. Betts C.B. Blucher A. Boniface C. Bucher E. Burlingame E. An omic and multidimensional spatial atlas from serial biopsies of an evolving metastatic breast cancer Cell Rep. Med. 3 2022 100525 10.1016/j.xcrm.2022.100525 35243422
10 Schurch C.M. Bhate S.S. Barlow G.L. Phillips D.J. Noti L. Zlobec I. Chu P. Black S. Demeter J. McIlwain D.R. Coordinated Cellular Neighborhoods Orchestrate Antitumoral Immunity at the Colorectal Cancer Invasive Front Cell 182 2020 1341 1359 10.1016/j.cell.2020.07.005 32763154
11 Gurel B. Lucia M.S. Thompson I.M. Jr. Goodman P.J. Tangen C.M. Kristal A.R. Parnes H.L. Hoque A. Lippman S.M. Sutcliffe S. Chronic inflammation in benign prostate tissue is associated with high-grade prostate cancer in the placebo arm of the prostate cancer prevention trial Cancer Epidemiol. Biomarkers Prev. 23 2014 847 856 10.1158/1055-9965.EPI-13-1126 24748218
12 De Marzo A.M. Platz E.A. Sutcliffe S. Xu J. Grönberg H. Drake C.G. Nakai Y. Isaacs W.B. Nelson W.G. Inflammation in prostate carcinogenesis Nat. Rev. Cancer 7 2007 256 269 10.1038/nrc2090 17384581
13 Bhandari V. Hoey C. Liu L.Y. Lalonde E. Ray J. Livingstone J. Lesurf R. Shiah Y.J. Vujcic T. Huang X. Molecular landmarks of tumor hypoxia across cancer types Nat. Genet. 51 2019 308 318 10.1038/s41588-018-0318-2 30643250
14 Barron D.A. Rowley D.R. The reactive stroma microenvironment and prostate cancer progression Endocr. Relat. Cancer 19 2012 R187 R204 10.1530/ERC-12-0085 22930558
15 Cioni B. Zwart W. Bergman A.M. Androgen receptor moonlighting in the prostate cancer microenvironment Endocr. Relat. Cancer 25 2018 R331 R349 10.1530/ERC-18-0042 29618577
16 Leach D.A. Buchanan G. Stromal Androgen Receptor in Prostate Cancer Development and Progression Cancers 9 2017 10 10.3390/cancers9010010
17 Eksi S.E. Chitsazan A. Sayar Z. Thomas G.V. Fields A.J. Kopp R.P. Spellman P.T. Adey A.C. Epigenetic loss of heterogeneity from low to high grade localized prostate tumours Nat. Commun. 12 2021 7292 10.1038/s41467-021-27615-8 34911933
18 Lin J.R. Wang S. Coy S. Chen Y.A. Yapp C. Tyler M. Nariya M.K. Heiser C.N. Lau K.S. Santagata S. Sorger P.K. Multiplexed 3D atlas of state transitions and immune interaction in colorectal cancer Cell 186 2023 363 381 10.1016/j.cell.2022.12.028 36669472
19 Meurs P. Galvin R. Fanning D.M. Fahey T. Prognostic value of the CAPRA clinical prediction rule: a systematic review and meta-analysis BJU Int. 111 2013 427 436 10.1111/j.1464-410X.2012.11400.x 22882877
20 Stringer C. Wang T. Michaelos M. Pachitariu M. Cellpose: a generalist algorithm for cellular segmentation Nat. Methods 18 2021 100 106 10.1038/s41592-020-01018-x 33318659
21 Burlingame E.A. Eng J. Thibault G. Chin K. Gray J.W. Chang Y.H. Toward reproducible, scalable, and robust data analysis across multiplex tissue imaging platforms Cell Rep. Methods 1 2021 100053 10.1016/j.crmeth.2021.100053 34485971
22 Becht E. McInnes L. Healy J. Dutertre C.A. Kwok I.W.H. Ng L.G. Ginhoux F. Newell E.W. Dimensionality reduction for visualizing single-cell data using UMAP. Nat Biotechnol Nature Biotechnol. 37 2019 38 44 10.1038/nbt.4314
23 Alinezhad S. Väänänen R.M. Ochoa N.T. Vertosick E.A. Bjartell A. Boström P.J. Taimen P. Pettersson K. Global expression of AMACR transcripts predicts risk for prostate cancer - a systematic comparison of AMACR protein and mRNA expression in cancerous and noncancerous prostate BMC Urol. 16 2016 1 10 10.1186/s12894-016-0128-8 26762220
24 Ananthanarayanan V. Deaton R.J. Yang X.J. Pins M.R. Gann P.H. Alpha-methylacyl-CoA racemase (AMACR) expression in normal prostatic glands and high-grade prostatic intraepithelial neoplasia (HGPIN): association with diagnosis of prostate cancer Prostate 63 2005 341 346 10.1002/pros.20196 15602744
25 Houlahan K.E. Shiah Y.J. Gusev A. Yuan J. Ahmed M. Shetty A. Ramanand S.G. Yao C.Q. Bell C. O'Connor E. Genome-wide germline correlates of the epigenetic landscape of prostate cancer Nat. Med. 25 2019 1615 1626 10.1038/s41591-019-0579-z 31591588
26 Sinha A. Huang V. Livingstone J. Wang J. Fox N.S. Kurganovs N. Ignatchenko V. Fritsch K. Donmez N. Heisler L.E. The Proteogenomic Landscape of Curable Prostate Cancer Cancer Cell 35 2019 414 427.e6 10.1016/j.ccell.2019.02.005 30889379
27 WaldhÖR T. The Spatial Autocorrelation Coefficient Moran'si under Heteroscedasticity Stat. Med. 15 1996 887 892 10.1002/(sici)1097-0258(19960415)15 8861157
28 Johansson A. Rudolfsson S. Hammarsten P. Halin S. Pietras K. Jones J. Stattin P. Egevad L. Granfors T. Wikström P. Bergh A. Mast cells are novel independent prognostic markers in prostate cancer and represent a target for therapy Am. J. Pathol. 177 2010 1031 1041 10.2353/ajpath.2010.100070 20616342
29 Fukui M. Whittlesey K. Metcalfe D.D. Dastych J. Human mast cells express the hyaluronic-acid-binding isoform of CD44 and adhere to hyaluronic acid Clin. Immunol. 94 2000 173 178 10.1006/clim.1999.4830 10692236
30 Wang W. Zou R. Qiu Y. Liu J. Xin Y. He T. Qiu Z. Interaction Networks Converging on Immunosuppressive Roles of Granzyme B: Special Niches Within the Tumor Microenvironment Front. Immunol. 12 2021 670324 10.3389/fimmu.2021.670324 33868318
31 Blatner N.R. Bonertz A. Beckhove P. Cheon E.C. Krantz S.B. Strouch M. Weitz J. Koch M. Halverson A.L. Bentrem D.J. Khazaie K. In colorectal cancer mast cells contribute to systemic regulatory T-cell dysfunction Proc. Natl. Acad. Sci. USA 107 2010 6430 6435 10.1073/pnas.0913683107 20308560
32 Saadalla A. Lima M.M. Tsai F. Osman A. Singh M.P. Linden D.R. Dennis K.L. Haeryfar S.M.M. Gurish M.F. Gounari F. Khazaie K. Cell Intrinsic Deregulated ss-Catenin Signaling Promotes Expansion of Bone Marrow Derived Connective Tissue Type Mast Cells, Systemic Inflammation, and Colon Cancer Front. Immunol. 10 2019 2777 10.3389/fimmu.2019.02777
33 Palla G. Spitzer H. Klein M. Fischer D. Schaar A.C. Kuemmerle L.B. Rybakov S. Ibarra I.L. Holmberg O. Virshup I. Squidpy: a scalable framework for spatial omics analysis Nat. Methods 19 2022 171 178 10.1038/s41592-021-01358-2 35102346
34 Ridge S.M. Sullivan F.J. Glynn S.A. Mesenchymal stem cells: key players in cancer progression Mol. Cancer 16 2017 31 10.1186/s12943-017-0597-8
35 Jung Y. Kim J.K. Shiozawa Y. Wang J. Mishra A. Joseph J. Berry J.E. McGee S. Lee E. Sun H. Recruitment of mesenchymal stem cells into prostate tumours promotes metastasis Nat. Commun. 4 2013 1795 10.1038/ncomms2766 23653207
36 Lai K.P. Yamashita S. Huang C.K. Yeh S. Chang C. Loss of stromal androgen receptor leads to suppressed prostate tumourigenesis via modulation of pro-inflammatory cytokines/chemokines EMBO Mol. Med. 4 2012 791 807 10.1002/emmm.201101140 22745041
37 Pederzoli F. Raffo M. Pakula H. Ravera F. Nuzzo P.V. Loda M. Stromal cells in prostate cancer pathobiology: friends or foes? Br. J. Cancer 128 2023 930 939 10.1038/s41416-022-02085-x 36482187
38 Tang Q. Cheng B. Dai R. Wang R. The Role of Androgen Receptor in Cross Talk Between Stromal Cells and Prostate Cancer Epithelial Cells Front. Cell Dev. Biol. 9 2021 729498 10.3389/fcell.2021.729498 34692685
39 Vidotto T. Saggioro F.P. Jamaspishvili T. Chesca D.L. Picanço de Albuquerque C.G. Reis R.B. Graham C.H. Berman D.M. Siemens D.R. Squire J.A. Koti M. PTEN-deficient prostate cancer is associated with an immunosuppressive tumor microenvironment mediated by increased expression of IDO1 and infiltrating FoxP3+ T regulatory cells Prostate 79 2019 969 979 10.1002/pros.23808 30999388
40 Bravo Gonzalez-Blas C. Minnoye L. Papasokrati D. Aibar S. Hulselmans G. Christiaens V. Davie K. Wouters J. Aerts S. cisTopic: cis-regulatory topic modeling on single-cell ATAC-seq data Nat. Methods 16 2019 397 400 10.1038/s41592-019-0367-1 30962623
41 Fang R. Preissl S. Li Y. Hou X. Lucero J. Wang X. Motamedi A. Shiau A.K. Zhou X. Xie F. Comprehensive analysis of single cell ATAC-seq data with SnapATAC Nat. Commun. 12 2021 1337 10.1038/s41467-021-21583-9 33637727
42 McLean C.Y. Bristor D. Hiller M. Clarke S.L. Schaar B.T. Lowe C.B. Wenger A.M. Bejerano G. GREAT improves functional interpretation of cis-regulatory regions Nat. Biotechnol. 28 2010 495 501 10.1038/nbt.1630 20436461
43 Zheng R. Wan C. Mei S. Qin Q. Wu Q. Sun H. Chen C.H. Brown M. Zhang X. Meyer C.A. Liu X.S. Cistrome Data Browser: expanded datasets and new tools for gene regulatory analysis Nucleic Acids Res. 47 2019 D729 D735 10.1093/nar/gky1094 30462313
44 Zhang Y. Guan X.Y. Jiang P. Cytokine and Chemokine Signals of T-Cell Exclusion in Tumors Front. Immunol. 11 2020 594609 10.3389/fimmu.2020.594609 33381115
45 Togashi Y. Shitara K. Nishikawa H. Regulatory T cells in cancer immunosuppression - implications for anticancer therapy Nat. Rev. Clin. Oncol. 16 2019 356 371 10.1038/s41571-019-0175-7 30705439
46 Tan H.L. Sood A. Rahimi H.A. Wang W. Gupta N. Hicks J. Mosier S. Gocke C.D. Epstein J.I. Netto G.J. Rb loss is characteristic of prostatic small cell neuroendocrine carcinoma Clin. Cancer Res. 20 2014 890 903 10.1158/1078-0432.CCR-13-1982 24323898
47 Becerra-Diaz M. Strickland A.B. Keselman A. Heller N.M. Androgen and Androgen Receptor as Enhancers of M2 Macrophage Polarization in Allergic Lung Inflammation J. Immunol. 201 2018 2923 2933 10.4049/jimmunol.1800352 30305328
48 Stultz J. Fong L. How to turn up the heat on the cold immune microenvironment of metastatic prostate cancer Prostate Cancer Prostatic Dis. 24 2021 697 717 10.1038/s41391-021-00340-5 33820953
49 Sfanos K.S. Immune cell infiltrates and prognosis in localized prostate cancer(dagger) J. Pathol. 256 2022 135 138 10.1002/path.5817 34694014
50 Fleischmann A. Schlomm T. Köllermann J. Sekulic N. Huland H. Mirlacher M. Sauter G. Simon R. Erbersdobler A. Immunological microenvironment in prostate cancer: high mast cell densities are associated with favorable tumor characteristics and good prognosis Prostate 69 2009 976 981 10.1002/pros.20948 19274666
51 Flammiger A. Bayer F. Cirugeda-Kühnert A. Huland H. Tennstedt P. Simon R. Minner S. Bokemeyer C. Sauter G. Schlomm T. Trepel M. Intratumoral T but not B lymphocytes are related to clinical outcome in prostate cancer APMIS 120 2012 901 908 10.1111/j.1600-0463.2012.02924.x 23009114
52 Zhao S.G. Lehrer J. Chang S.L. Das R. Erho N. Liu Y. Sjöström M. Den R.B. Freedland S.J. Klein E.A. The Immune Landscape of Prostate Cancer and Nomination of PD-L2 as a Potential Therapeutic Target J. Natl. Cancer Inst. 111 2019 301 310 10.1093/jnci/djy141 30321406
53 Andersen L.B. Nørgaard M. Rasmussen M. Fredsøe J. Borre M. Ulhøi B.P. Sørensen K.D. Immune cell analyses of the tumor microenvironment in prostate cancer highlight infiltrating regulatory T cells and macrophages as adverse prognostic factors J. Pathol. 255 2021 155 165 10.1002/path.5757 34255349
54 Kaur H.B. Guedes L.B. Lu J. Maldonado L. Reitz L. Barber J.R. De Marzo A.M. Tosoian J.J. Tomlins S.A. Schaeffer E.M. Association of tumor-infiltrating T-cell density with molecular subtype, racial ancestry and clinical outcomes in prostate cancer Mod. Pathol. 31 2018 1539 1552 10.1038/s41379-018-0083-x 29849114
55 Davidsson S. Ohlson A.L. Andersson S.O. Fall K. Meisner A. Fiorentino M. Andrén O. Rider J.R. CD4 helper T cells, CD8 cytotoxic T cells, and FOXP3(+) regulatory T cells with respect to lethal prostate cancer Mod. Pathol. 26 2013 448 455 10.1038/modpathol.2012.164 23041830
56 Hempel Sullivan H. Heaphy C.M. Kulac I. Cuka N. Lu J. Barber J.R. De Marzo A.M. Lotan T.L. Joshu C.E. Sfanos K.S. High Extratumoral Mast Cell Counts Are Associated with a Higher Risk of Adverse Prostate Cancer Outcomes Cancer Epidemiol. Biomarkers Prev. 29 2020 668 675 10.1158/1055-9965.EPI-19-0962 31932412
57 Hempel H.A. Cuka N.S. Kulac I. Barber J.R. Cornish T.C. Platz E.A. De Marzo A.M. Sfanos K.S. Low Intratumoral Mast Cells Are Associated With a Higher Risk of Prostate Cancer Recurrence Prostate 77 2017 412 424 10.1002/pros.23280 27868214
58 Khazaie K. Blatner N.R. Khan M.W. Gounari F. Gounaris E. Dennis K. Bonertz A. Tsai F.N. Strouch M.J. Cheon E. The significant role of mast cells in cancer Cancer Metastasis Rev. 30 2011 45 60 10.1007/s10555-011-9286-z 21287360
59 Valdman A. Jaraj S.J. Compérat E. Charlotte F. Roupret M. Pisa P. Egevad L. Distribution of Foxp3-CD4- and CD8-positive lymphocytic cells in benign and malignant prostate tissue APMIS 118 2010 360 365 10.1111/j.1600-0463.2010.02604.x 20477811
60 Hussein M.R.A. Al-Assiri M. Musalam A.O. Phenotypic characterization of the infiltrating immune cells in normal prostate, benign nodular prostatic hyperplasia and prostatic adenocarcinoma Exp. Mol. Pathol. 86 2009 108 113 10.1016/j.yexmp.2008.11.010 19111537
61 Fujii T. Shimada K. Asai O. Tanaka N. Fujimoto K. Hirao K. Konishi N. Immunohistochemical analysis of inflammatory cells in benign and precancerous lesions and carcinoma of the prostate Pathobiology 80 2013 119 126 10.1159/000342396 23328608
62 Nardone V. Botta C. Caraglia M. Martino E.C. Ambrosio M.R. Carfagno T. Tini P. Semeraro L. Misso G. Grimaldi A. Tumor infiltrating T lymphocytes expressing FoxP3, CCR7 or PD-1 predict the outcome of prostate cancer patients subjected to salvage radiotherapy after biochemical relapse Cancer Biol. Ther. 17 2016 1213 1220 10.1080/15384047.2016.1235666 27791459
63 Watanabe M. Kanao K. Suzuki S. Muramatsu H. Morinaga S. Kajikawa K. Kobayashi I. Nishikawa G. Kato Y. Zennami K. Increased infiltration of CCR4-positive regulatory T cells in prostate cancer tissue is associated with a poor prognosis Prostate 79 2019 1658 1665 10.1002/pros.23890 31390096
64 Lanciotti M. Masieri L. Raspollini M.R. Minervini A. Mari A. Comito G. Giannoni E. Carini M. Chiarugi P. Serni S. The role of M1 and M2 macrophages in prostate cancer in relation to extracapsular tumor extension and biochemical recurrence after radical prostatectomy BioMed Res. Int. 2014 2014 486798 10.1155/2014/486798 24738060
65 Erlandsson A. Carlsson J. Lundholm M. Fält A. Andersson S.O. Andrén O. Davidsson S. M2 macrophages and regulatory T cells in lethal prostate cancer Prostate 79 2019 363 369 10.1002/pros.23742 30500076
66 Sadasivan S.M. Chen Y. Gupta N.S. Han X. Bobbitt K.R. Chitale D.A. Williamson S.R. Rundle A.G. Tang D. Rybicki B.A. The interplay of growth differentiation factor 15 (GDF15) expression and M2 macrophages during prostate carcinogenesis Carcinogenesis 41 2020 1074 1082 10.1093/carcin/bgaa065 32614434
67 Baker E.A.G. Schapiro D. Dumitrascu B. Vickovic S. Regev A. In silico tissue generation and power analysis for spatial omics Nat. Methods 20 2023 424 431 10.1038/s41592-023-01766-6 36864197
68 Comito G. Giannoni E. Segura C.P. Barcellos-de-Souza P. Raspollini M.R. Baroni G. Lanciotti M. Serni S. Chiarugi P. Cancer-associated fibroblasts and M2-polarized macrophages synergize during prostate carcinoma progression Oncogene 33 2014 2423 2431 10.1038/onc.2013.191 23728338
69 Allavena P. Digifico E. Belgiovine C. Macrophages and cancer stem cells: a malevolent alliance Mol. Med. 27 2021 121 10.1186/s10020-021-00383-3 34583655
70 Zhao Y.B. Yang S.H. Shen J. Deng K. Li Q. Wang Y. Cui W. Ye H. Interaction between regulatory T cells and mast cells via IL-9 and TGF-beta production Oncol. Lett. 20 2020 360 10.3892/ol.2020.12224 33133260
71 Blatner N.R. Mulcahy M.F. Dennis K.L. Scholtens D. Bentrem D.J. Phillips J.D. Ham S. Sandall B.P. Khan M.W. Mahvi D.M. Expression of RORgt Marks a Pathogenic T Cell Subset in Human Colon Cancer Sci. Transl. Med. 4 2012 164ra159
72 Hempel Sullivan H. Maynard J.P. Heaphy C.M. Lu J. De Marzo A.M. Lotan T.L. Joshu C.E. Sfanos K.S. Differential mast cell phenotypes in benign versus cancer tissues and prostate cancer oncologic outcomes J. Pathol. 253 2021 415 426 10.1002/path.5606 33338262
73 Takano H. Nakazawa S. Shirata N. Tamba S. Furuta K. Tsuchiya S. Morimoto K. Itano N. Irie A. Ichikawa A. Involvement of CD44 in mast cell proliferation during terminal differentiation Lab. Invest. 89 2009 446 455 10.1038/labinvest.2008.159 19204665
74 Mueller N. Wicklein D. Eisenwort G. Jawhar M. Berger D. Stefanzl G. Greiner G. Boehm A. Kornauth C. Muellauer L. CD44 is a RAS/STAT5-regulated invasion receptor that triggers disease expansion in advanced mastocytosis Blood 132 2018 1936 1950 30018080
75 Guan X. Polesso F. Wang C. Sehrawat A. Hawkins R.M. Murray S.E. Thomas G.V. Caruso B. Thompson R.F. Wood M.A. Androgen receptor activity in T cells limits checkpoint blockade efficacy Nature 606 2022 791 796 10.1038/s41586-022-04522-6 35322234
76 Chen W. Beck I. Schober W. Brockow K. Effner R. Buters J.T.M. Behrendt H. Ring J. Human mast cells express androgen receptors but treatment with testosterone exerts no influence on IgE-independent mast cell degranulation elicited by neuromuscular blocking agents Exp. Dermatol. 19 2010 302 304 10.1111/j.1600-0625.2009.00969.x 19758318
77 Jia J. Zhang H. Zhang H. Liu W. Du H. Shu M. He L. AR facilitates YAP-TEAD interaction with the AM promoter to enhance mast cell infiltration into cutaneous neurofibroma Sci. Rep. 9 2019 19346 10.1038/s41598-019-56022-9 31852972
78 Dang Q. Li L. Xie H. He D. Chen J. Song W. Chang L.S. Chang H.C. Yeh S. Chang C. Anti-androgen enzalutamide enhances prostate cancer neuroendocrine (NE) differentiation via altering the infiltrated mast cells--> androgen receptor (AR)--> miRNA32 signals Mol. Oncol. 9 2015 1241 1251 10.1016/j.molonc.2015.02.010 25817444
79 Kron K.J. Murison A. Zhou S. Huang V. Yamaguchi T.N. Shiah Y.J. Fraser M. van der Kwast T. Boutros P.C. Bristow R.G. Lupien M. TMPRSS2-ERG fusion co-opts master transcription factors and activates NOTCH signaling in primary prostate cancer Nat. Genet. 49 2017 1336 1345 10.1038/ng.3930 28783165
80 Rajaram S. Heinrich L.E. Gordan J.D. Avva J. Bonness K.M. Witkiewicz A.K. Malter J.S. Atreya C.E. Warren R.S. Wu L.F. Altschuler S.J. Sampling strategies to capture single-cell heterogeneity Nat. Methods 14 2017 967 970 10.1038/nmeth.4427 28869755
81 Liu Y. Wang J. Horton C. Yu C. Knudsen B. Stefanson J. Hu K. Stefanson O. Green J. Guo C. Stromal AR inhibits prostate tumor progression by restraining secretory luminal epithelial cells Cell Rep. 39 2022 110848 10.1016/j.celrep.2022.110848 35613593
82 Malinen M. Niskanen E.A. Kaikkonen M.U. Palvimo J.J. Crosstalk between androgen and pro-inflammatory signaling remodels androgen receptor and NF-kappaB cistrome to reprogram the prostate cancer cell transcriptome Nucleic Acids Res. 45 2017 619 630 10.1093/nar/gkw855 27672034
83 Zhang L. Altuwaijri S. Deng F. Chen L. Lal P. Bhanot U.K. Korets R. Wenske S. Lilja H.G. Chang C. NF-kappaB regulates androgen receptor expression and prostate cancer growth Am. J. Pathol. 175 2009 489 499 10.2353/ajpath.2009.080727 19628766
84 Nirmal A.J. Maliga Z. Vallius T. Quattrochi B. Chen A.A. Jacobson C.A. Pelletier R.J. Yapp C. Arias-Camison R. Chen Y.A. The Spatial Landscape of Progression and Immunoediting in Primary Melanoma at Single-Cell Resolution Cancer Discov. 12 2022 1518 1541 10.1158/2159-8290.CD-21-1357 35404441
85 Pourmaleki M. Jones C.J. Ariyan C.E. Zeng Z. Pirun M. Navarrete D.A. Li Y. Zhang M. Nandakumar S. Campos C. Tumor MHC Class I Expression Associates with Intralesional IL2 Response in Melanoma Cancer Immunol. Res. 10 2022 303 313 10.1158/2326-6066.CIR-21-1083 35013003
86 Song H. Bucher S. Rosenberg K. Tsui M. Burhan D. Hoffman D. Cho S.J. Rangaswami A. Breese M. Leung S. Single-cell analysis of hepatoblastoma identifies tumor signatures that predict chemotherapy susceptibility using patient-specific tumor spheroids Nat. Commun. 13 2022 4878 10.1038/s41467-022-32473-z 36008377
87 Cooperberg M.R. Pasta D.J. Elkin E.P. Litwin M.S. Latini D.M. Du Chane J. Carroll P.R. San Francisco Cancer of the Prostate Risk Assessment score: a straightforward and reliable preoperative predictor of disease recurrence after radical prostatectomy J. Urol. 173 2005 1938 1942 10.1097/01.ju.0000158155.33890.e7 15879786
88 Lin J.R. Izar B. Wang S. Yapp C. Mei S. Shah P.M. Santagata S. Sorger P.K. Highly multiplexed immunofluorescence imaging of human tissues and tumors using t-CyCIF and conventional optical microscopes Elife 7 2018 e31657 29993362
89 Eng J. Bucher E. Hu Z. Zheng T. Gibbs S.L. Chin K. Gray J.W. A framework for multiplex imaging optimization and reproducible analysis Commun. Biol. 5 2022 438 10.1038/s42003-022-03368-y 35545666
90 Adalsteinsson V.A. Ha G. Freeman S.S. Choudhury A.D. Stover D.G. Parsons H.A. Gydush G. Reed S.C. Rotem D. Rhoades J. Scalable whole-exome sequencing of cell-free DNA reveals high concordance with metastatic tumors Nat. Commun. 8 2017 1324 10.1038/s41467-017-00965-y 29109393
91 Cusanovich D.A. Daza R. Adey A. Pliner H.A. Christiansen L. Gunderson K.L. Steemers F.J. Trapnell C. Shendure J. Multiplex single-cell profiling of chromatin accessibility by combinatorial cellular indexing Science 348 2015 910 914 10.1126/science.aab1601 25953818
92 Buenrostro J.D. Wu B. Litzenburger U.M. Ruff D. Gonzales M.L. Snyder M.P. Chang H.Y. Greenleaf W.J. Single-cell chromatin accessibility reveals principles of regulatory variation Nature 523 2015 486 490 10.1038/nature14590 26083756
93 Chang Y.H. Chin K. Thibault G. Eng J. Burlingame E. Gray J.W. RESTORE: Robust intEnSiTy nORmalization mEthod for multiplexed imaging Commun. Biol. 3 2020 111 10.1038/s42003-020-0828-1 32152447
94 Jackson H.W. Fischer J.R. Zanotelli V.R.T. Ali H.R. Mechera R. Soysal S.D. Moch H. Muenst S. Varga Z. Weber W.P. Bodenmiller B. The single-cell pathology landscape of breast cancer Nature 578 2020 615 620 10.1038/s41586-019-1876-x 31959985
95 Vallat R. Pingouin: statistics in Python J. Open Source Softw. 3 2018 1026 10.21105/joss.01026
96 Goltsev Y. Samusik N. Kennedy-Darling J. Bhate S. Hale M. Vazquez G. Black S. Nolan G.P. Deep Profiling of Mouse Splenic Architecture with CODEX Multiplexed Imaging Cell 174 2018 968 981.e15 10.1016/j.cell.2018.07.010 30078711
