Computational Mapping of Reactive Stroma in Prostate Cancer Yields Interpretable, Prognostic Biomarkers

TL;DR

PROTAS addresses a gap in prostate cancer risk stratification by quantifying reactive stroma (RS) from routine H&E slides and linking RS morphology to underlying biology and prognosis. It uses a frozen foundation encoder (UNI) plus a two-layer MLP, with severity-aware soft labels and entropy regularization, and extends to needle biopsies through domain-adversarial training. Across internal and external cohorts, RS features improve biochemical recurrence prediction beyond standard baselines and show higher reproducibility than expert pathologists. Transcriptomic and spatial analyses corroborate a contractile ECM remodeling RS program, supporting PROTAS as an interpretable, scalable biomarker that complements glandular-centric grading and informs precision management.

Abstract

Current histopathological grading of prostate cancer relies primarily on glandular architecture, largely overlooking the tumor microenvironment. Here, we present PROTAS, a deep learning framework that quantifies reactive stroma (RS) in routine hematoxylin and eosin (H&E) slides and links stromal morphology to underlying biology. PROTAS-defined RS is characterized by nuclear enlargement, collagen disorganization, and transcriptomic enrichment of contractile pathways. PROTAS detects RS robustly in the external Prostate, Lung, Colorectal, and Ovarian (PLCO) dataset and, using domain-adversarial training, generalizes to diagnostic biopsies. In head-to-head comparisons, PROTAS outperforms pathologists for RS detection, and spatial RS features predict biochemical recurrence independently of established prognostic variables (c-index 0.80). By capturing subtle stromal phenotypes associated with tumor progression, PROTAS provides an interpretable, scalable biomarker to refine risk stratification.

Figures (17)

• Figure 1: Study Pipeline (a). Cohort curation and slide selection. For each RP case, pathology reports were parsed to identify the dominant tumor focus and location. An automated cancer-prediction model produced tumor masks, which were cross-checked against the report to select one representative whole-mount H&E WSI per patient, resulting in a final cohort of $N = 255$ (b). Stroma-patch isolation and labeling. Gland segmentation (internal model) and cell semantic segmentation (CellViT) were used to remove glandular epithelium and isolate stroma-only patches. Slides were tiled at 40$\times$ into 256$\times$256 pixel patches, then stain-normalized. Positive stroma patches were sampled from within the dominant tumor region (report-concordant), whereas negative stroma patches were sampled from $\geq$ 10 mm away from any tumor focus. (c). Model training and inference. Patients were split at the patient level (train/validation/test: 64.5%/11.1%/24.4%). A pretrained pathology feature encoder generated fixed embeddings for stroma patches; a two-layer neural network (with batch normalization and dropout) was trained with severity-aware soft labels and an entropy (confidence-penalty) term. At inference, Mont Carlo Dropout provided mean RS probability and epistemic uncertainty for each RS prediction.
• Figure 2: Slide-level RS features and spatial maps. (a). Distributions of six representative slide-level RS features comparing low-grade ($\leq$GG2) versus high-grade ($\geq$GG3): entropy of RS-positive probabilities, hotspot Laplacian mean, number of hotspot-hotspot connections, graph diameter (hops), 95th-percentile RS probability within the 1-2 mm ring outside tumor, and maximum RS rate across 1-mm distance bins. Group differences were assessed with Wilcoxon rank-sum tests with Bonferroni adjustment. (b). The same six features stratified by PI-RADS score (2-4 vs. 5). (c). Heatmap of per-group medians for all slide-level features (left: $\leq$GG2 vs. $\geq$GG3; right: PI-RADS 2-4 vs. 5); deeper magenta indicates a higher median value. (d). Example whole-mount probability maps over stroma-only patches showing PROTAS predictions; warmer colors denote higher RS probability (scale 0-1). Yellow outlines indicate tumor regions.
• Figure 3: Prognostic Modeling with RS features, tumor features, and UCSF-CAPRA (a). Model comparison across four metrics. Bars show point estimates with 95% bootstrap CIs (5,000 replicates) for concordance index and time-dependent area under receiver operating characteristic (AUROC) at 12, 24, and 60 months. Asterisks denote significant pairwise differences (bootstrap $p<0.05$) across models. (b). Multivariable Cox model including all feature sets (RS features, tumor features, UCSF-CAPRA): hazard ratios with 95% CIs on a logarithmic scale; the dashed line indicates HR = 1. Feature definitions and binning/standardization follow Methods. (c). Kaplan-Meier curves stratified by quartiles of each model's partial hazard scores; Adjusted Log-rank $p$-values are reported on each panel.
• Figure 4: Handcrafted feature models and interpretability. (a). Test set performance of four classifiers (Random Forest, XGBoost, naive Bayes, two-layer neural network) trained on nuclear-only, collagen-only, and combined feature sets. Bars show AUROC, accuracy, specificity, and sensitivity; colors denote model type and hatching denotes feature set. (b). Odds ratios (95% CIs; log-scale) from a logistic GLM with slide fixed effects for nuclear and collagen features that remained significant after Bonferroni correction; the dashed line marks OR = 1. (c). Top-15 features ranked by mean absolute SHAP value for the MLP trained on combined features. Points show per-patch SHAP values (impact on model); color encodes the feature value (pink = high, blue = low). (d). SHAP waterfall plots for representative RS-positive and RS-negative patches, illustrating how individual feature values shift the prediction from model baseline $E[f(X)]$ to the final output (positive contributions increase RS probability).
• Figure 5: TCGA-PRAD differential gene expression (a) Volcano plots (left) RS-rich vs RS-lower patients (RS-rich defined as the top quartile of the patient-level RS score: proportion of positive patches $\geq$ 75th percentile) and (right) low-grade ($\leq$GG2) vs. high-grade ($\geq$GG3) patients. Points show $\text{log}_2$ fold change (x-axis) versus $-\text{log}_{10}$ adjusted $p$-value (y-axis; DESeq2 with the Benjamini-Hochberg method). Vertical dashed lines mark $|\text{log}_2 \text{FC}| = 1$; the horizontal dashed line marks FDR = $0.05$.