This retrospective multicenter cohort research included patients experiencing radical resection for CRC, drawing from three independent cohorts: TCGA-COAD, TCGA-READ, and real-world (HMUCH) cohort. The TCGA-COAD and TCGA-READ cohorts were combined to form a unified meta-cohort, designated as TCGA-CRC. The training cohort consisted of 485 consecutive patients treated at HMUCH between January 2012 and December 2013.

Patients were eligible for inclusion if they met the following criteria: (1) histopathologically confirmed CRC with R0 surgical margins; (2) survival of at least 90 days postoperatively to minimize bias from surgical quality (16, 17); (3) no prior history of malignancy; and (4) availability of complete clinical, pathological, and follow-up records. Individuals who received neoadjuvant therapy were excluded, as such treatment may alter tissue morphology in H&E-stained slides and affect the prognostic assessment.

A total of 398 CRC patients meeting the same eligibility criteria were obtained from TCGA database through the National Cancer Institute’s Genomic Data Commons (https://gdc.cancer.gov/). These cases, which included complete prognostic data and high-quality digital H&E-stained histopathological images, served as the validation cohort.

Baseline clinical and pathological characteristics were comprehensively collected, including patient age, sex, tumor site, invasion depth, perineural invasion, lymphovascular invasion, vascular invasion, lymph node involvement, TNM stage, follow-up information (duration and survival status), and receipt of postoperative adjuvant chemotherapy.

The determination of follow-up duration was performed from the surgery date to the most recent follow-up, with survival status documented at the final visit. Overall survival (OS) was set as the interval between either last follow-up or death and surgery. Disease-free survival (DFS) refers to the time from surgery to the first occurrence of recurrence at any site or death from any cause, whichever occurs earlier.

Slides from the HMUCH-CRC cohort were prepared through routine histopathological processing involving fixation in 4% neutral formaldehyde, paraffin embedding, 4 μm sectioning, and H&E staining. TNM staging was subsequently reassessed according to the 8th edition criteria of the American Joint Committee on Cancer (AJCC). For each case, representative sections illustrating the invasion depth were carefully selected. Following quality control, the slides were scanned using an Aperio AT2 scanner (Leica Biosystems, Germany) at 20× optical magnification (0.5 μm/pixel). The resulting digital images were stored in SVS format and managed using Aperio ImageScope software (version 12.4.6).

To facilitate the processing of WSIs approaching 10 gigapixels in size, we first applied the OTSU thresholding algorithm to remove white background regions (18). Subsequently, we partitioned the non-background region into non-overlapping image patches measuring 512 × 512 pixels at a 20-fold optical magnification and recorded their respective locations, resulting in over 7.7 million patches. Note that the batch size is 32, the initial learning rate is 0.01, the cosine decay optimizer is SGD and the momentum is 0.9. Additionally, we applied the Macenko method to normalize the color of small tiles (19), followed by z-score normalization on RGB channels to achieve a standard normal distribution of image intensities.

To enhance model generalization, various data augmentation strategies, such as flipping, mirroring, blurring, mild color perturbations, and progressive sprinkling, were randomly applied to the images in the validation and training sets. Notably, no augmentation was conducted on test images.

In this study, we designed an advanced deep learning framework to address the complexity and heterogeneity inherent in large-scale tumor histopathology images. The model adopts a two-stage architecture, beginning with patch-level inference and subsequently integrating patch probabilities through a multi-instance learning (MIL)-based feature fusion algorithm to generate WSI-level predictions.

During training, each image patch was assigned the same label corresponding to the patient’s 5-year survival status. For patch-level classification, we employed ResNet-18, an established convolutional neural network architecture renowned for its success in the ImageNet challenge, to estimate patch-level survival likelihoods. Model optimization was performed using softmax cross-entropy loss and mini-batch stochastic gradient descent (SGD).

To enhance generalizability across heterogeneous cohorts, we applied transfer learning by initializing model weights using pre-trained parameters from the ImageNet dataset. The learning rate was fine-tuned via a cosine annealing schedule, defined as

η m i n i = 0 indicates the minimum learning rate, and η m a x i = 0.01 represents the maximum learning rate. The term T i = 30 denotes the number of iteration epochs used in the model training. We also utilize transfer learning algorithms to ensure optimal model fitting, by fine-tuning the backbone component parameters when T c u r = 1 2 T i . The learning rate for the backbone component is defined as follows:

Following the model training, each patch was assigned a prediction label along with its corresponding probability. These patch-level likelihoods were then aggregated using a classifier to generate the WSI-level outcomes. To facilitate this process, we developed two distinct MIL pipelines: the Patch Likelihood Histogram (PALHI) method and the Bag of Words (BoW) method, inspired by histogram-based and vocabulary-based strategies, respectively.

In the PALHI pipeline, a histogram-based representation is used to quantify the distribution of patch-level likelihoods within each WSI. In contrast, the BoW approach encodes each patch as a floating-point value using term frequency–inverse document frequency (TF-IDF) with a feature vector representing the entire slide.

Using these two distinct pipelines, patch-level outputs were effectively transformed into WSI-level features. Each method contributed 101 probabilistic features and 2 categorical label features. These individual feature sets were then integrated through early fusion, resulting in a unified feature vector of 206 dimensions for subsequent analysis.

The least absolute shrinkage and selection operator (LASSO) Cox regression model that incorporates an L1 penalty for the toward-zero reduction of feature coefficients is a well-established method for survival analysis in high-dimensional settings (20, 21). Referred to as the tuning constant, the penalty parameter λ governs the penalty strength.

In this study, the 10-fold cross-validation with the minimum criteria was applied for obtaining the optimal λ by minimizing the partial likelihood deviance within the training set. This approach enabled the identification of key prognostic features and construction of a formula for calculating the pathomics signature. The application of the derived formula to the validation set helped compute the corresponding signature scores.

Through using the maximally selected rank statistics within the training cohort, the identification of the optimal cutoff for the pathomics signature was achieved, which was subsequently tested within the validation set. This threshold was adopted for the classification of the patients into high- and low-signature categories for prognostic evaluation. The differences in OS and DFS between both categories were analyzed using restricted mean survival time (RMST) metrics and Kaplan-Meier (K-M) survival curves (22).

The independent prognostic significance of the pathomics signature was assessed through performing univariate and multivariate Cox regression analyses. Subgroup heterogeneity was examined using interaction-based subgroup analysis. To evaluate the potential influence of unmeasured confounding, E-value analysis was conducted as a sensitivity assessment (23).

To measure discriminative performance, we calculated the concordance index (C-index) and the area under the receiver operating characteristic curve (AUROC). The agreement between predicted and observed survival probabilities was evaluated through applying calibration plots. The clinical utility of the pathomics model was further examined using decision curve analysis (DCA), which quantifies net benefit across varying decision thresholds (24).

To determine the added value of the pathomics signature beyond conventional TNM staging, we evaluated its impact on discrimination, calibration, clinical benefit, integrated discrimination improvement (IDI), net reclassification improvement (NRI), and prediction error curves (25).

To mitigate the interpretability limitations of deep learning models, we applied SHapley Additive exPlanations (SHAP), a method rooted in cooperative game theory, to quantify the contribution and relative importance of individual features to model outputs (26). This technique enables both global and instance-level interpretation of the predictions generated by the trained model.

The gradient-weighted class activation mapping (Grad-CAM) was utilized to produce heatmaps over selected image tiles for further exploration of prognostically relevant morphological patterns (27), highlighting crucial regions that influenced network predictions. This visualization technique utilized gradient information from the last convolutional layer of our deep learning network, providing a visual explanation that facilitates understanding and validating of the model’s decision-making process.

Transcriptomic profiles from TCGA cohort were retrieved with the TCGAbiolinks package (28). Gene set enrichment analysis (GSEA) was performed to infer the biological processes related to the pathomics signature (29). Additionally, pathway activity was quantified using gene set variation analysis (GSVA) via the GSVA package (30), allowing the identification of significantly enriched pathways across different patient subgroups. Functional interpretation was based on the well-curated “hallmark gene sets” (31).

Weighted correlation network analysis (WGCNA) was performed using the WGCNA package, which aims to identify the pathomics signature-related gene modules (32). The scale-free topology fitting index of 0.85 was set as the threshold to construct the signed weighted gene co-expression network. The minimum co-expression module size was set to 30, and the merge cut minimum module merge cut height was set to 0.25. A biweight midcorrelation coefficient (bicor) > 0.1 and P-value < 0.05 were selected as the thresholds to find gene modules significantly associated with the pathomics signature. Gene annotation enrichment analysis was performed using the clusterProfiler package (33).

Furthermore, according to the guideline for transcriptome-based cell-type quantification methods, we utilized the MCPcounter and xCell algorithms to quantify the proportions of specific immune and stromal cells within the CRC samples (34–36).

The comparison of the continuous variables with normal distributions were performed using unpaired two-sample t-tests, while the analysis of the non-normally distributed variables were achieved via the Mann-Whitney U test or Kruskal-Wallis test. The assessment of the categorical variables was achieved via either Fisher’s exact test or the chi-squared (χ²) test. Survival curves were generated through applying the K-M method and evaluated via the log-rank test. Univariate and multivariate associations were examined using Cox regression analysis with 95% confidence intervals (CIs) and hazard ratios (HRs). The associations between continuous variables were analyzed through the calculation of the Spearman rank correlation coefficients. All statistical analyses were conducted using R software (v4.0.5) and SPSS (v19.0). Deep learning experiments were implemented in Python (v3.7.12). All tests were two-sided, with P-value < 0.05 considered statistically significant.

Detailed clinicopathological features of patients from the training cohort (n = 485) and the validation cohort (n = 398) are summarized in Supplementary Table S1 . Across all 883 patients, the median age was 62 years with the interquartile range (IQR) of 54–71, where males accounted for 54.9% (485/883) of the population. The majority (86.2%, 761/883) was diagnosed at stage II or III. In the training set, the median follow-up period was 72.5 months (IQR: 56.47–121.23), with 5-year DFS and OS rates of 75.51% and 81.28%, respectively. In contrast, the validation cohort had a shorter median follow-up of 24.33 months (IQR: 15.24–36.53), with corresponding 5-year DFS and OS rates of 61.89% and 70.97%. The differences observed in the clinicopathological profiles between cohorts reflect real-world clinical diversity, thereby enhancing the generalizability of our results.

The development framework for the pathomics signature is depicted in Figure 1 . In the training cohort, a LASSO-Cox regression model with 10-fold cross-validation was employed to construct the signature. Using the optimal penalty parameter λ ( Supplementary Figure S1 ), eight selected pathomics features were integrated into a composite risk score. The final formula for calculating the pathomics signature is as follows:

The optimal cutoff point, identified based on the maximum standardized log-rank statistic, was 0.1139008. Patients in both the training and validation cohorts were stratified into high- and low-signature groups accordingly. Associations between the pathomics signature and clinicopathological characteristics are presented in Supplementary Table S2 . Notably, a potential correlation was observed between the signature and lymph node counts.

Supplementary Figure S2 illustrates the distribution of pathomics signature values by survival status along with selected feature profiles, indicating a positive association between elevated signature scores and a higher risk of recurrence or mortality. K-M survival analysis ( Figure 2 ) demonstrated significant differences in both OS and DFS between the low- and high-signature groups in the training and validation cohorts.

RMST analysis revealed a sustained survival advantage for patients with low pathomics signature scores across multiple time points, with the magnitude of the benefit increasing over time ( Table 1 ). Specifically, the low-signature group exhibited an OS advantage of approximately 2 months at year 3, 8 months at year 5, and a notable 15 months by year 7 when compared to the high-signature group.

Univariate and multivariate Cox regression analyses confirmed the pathomics signature as an independent predictor of both DFS and OS in the training cohort ( Table 2 ). Consistent findings were observed in the validation cohort ( Supplementary Table S3 ). To assess the robustness of these associations against potential unmeasured confounding, E-value sensitivity analyses were conducted based on adjusted HRs in both cohorts ( Supplementary Table S4 ).

Stratified analyses based on clinicopathological variables demonstrated that the pathomics signature remained a significant prognostic marker across all subgroups, except for patients with perineural invasion in the validation cohort ( Table 3 ). A potential interaction between age and lymph node harvest was suggested by the subgroup difference testing. No other significant interaction effects were observed, thereby supporting the overall robustness of the pathomics signature as a prognostic factor.

Time-dependent receiver operating characteristic (ROC) curves demonstrated that the pathomics signature achieved favorable predictive performance for 3-, 5-, and 7-year OS and DFS in both the training and validation cohorts ( Supplementary Figure S3 ). The corresponding calibration plots further confirmed a strong agreement between the predicted and observed survival probabilities across the same time intervals ( Supplementary Figure S4 ).

Moreover, decision curve analysis (DCA) showed that incorporating the pathomics signature into prognostic assessment yielded greater net clinical benefit than either the “treat-all” or “treat-none” strategies in both cohorts ( Supplementary Figure S5 ), supporting its potential for real-world clinical application.

The combined model, which were based on the combination of the pathomics signature and the TNM staging system, exhibited a significantly higher C-Index than the TNM stage, and these results could also be found in the validation cohort ( Supplementary Table S5 ).

Furthermore, the AUROCs of the 3 models also confirmed the superior discrimination ability of the combined models for estimating DFS and OS in the training and validation cohorts ( Supplementary Figure S6 ). Additionally, compared with the TNM stage models, the combined models were the most accurate models ( Supplementary Figure S7 ) and showed greater net benefits across most of the range of reasonable threshold probabilities ( Supplementary Figure S8 ).

Finally, the combined model showed a significant NRI and IDI for prognosis estimation compared with the TNM stage model ( Supplementary Table S6 ), indicating that the pathomics signature could provide additional prognostic value to the TNM staging system for CRC.

To evaluate the predictive utility of the pathomics signature in the context of adjuvant chemotherapy, we analyzed its association with survival outcomes in stage II and III CRC patients stratified by postoperative adjuvant chemotherapy status. In both the training and validation cohorts, adjuvant chemotherapy significantly improved OS and DFS in these subgroups ( Supplementary Figure S9 ). Furthermore, the pathomics signature demonstrated a significant correlation with OS and DFS, regardless of whether the patients received adjuvant therapy ( Supplementary Figure S10 ).

Among the patients in the low-pathomics signature group, adjuvant chemotherapy was significantly associated with improved OS and DFS. In contrast, this survival benefit was not observed in the high-pathomics signature group ( Figure 3 ). Further interaction analysis revealed a significant effect modification, indicating that individuals with low pathomics signature scores derived greater benefits from adjuvant chemotherapy than those with high scores ( Supplementary Table S7 ).

SHAP values were used to interpret the contribution of individual features to the model predictions. As illustrated in the SHAP summary plot ( Figure 4A ), HistogramBoWProb_0.15 emerged as the most influential feature, closely followed by HistogramBoWProb_0.66, HistogramBoWProb_0.72, and BoWProb_008. In contrast, BoWProb_068 contributed the least among the eight pathomics features.

To further interpret the spatial features linked to patient prognosis, prediction heatmaps were generated to highlight the key regions contributing to the model output ( Figure 4B ). High-risk cases are typically marked by dense tumor stroma, abundant tumor cells, and muscle tissue infiltration. In contrast, low-risk regions were predominantly characterized by normal mucosa, loose stroma, and inflammatory infiltration.

To further elucidate the model’s decision-making process, Grad-CAM was applied to extract informative visual cues. The top 500 most influential regions were selected to explore the dominant histopathological patterns associated with patient survival. By clustering the highest-ranked image patches, four distinct histological clusters were identified using a random tree algorithm ( Figure 4C ).

Subsequently, expert pathologists reviewed and annotated the representative regions identified using the model. Key histological components, including tumor cells, adipocytes, fibrous tissue, and stroma, were highlighted in red ( Figure 4D ). These features appeared to be closely associated with an elevated risk of recurrence and mortality, offering a morphological interpretation of the predictive elements underlying the pathomics signature.

GSEA was initially conducted to investigate potential biological mechanisms associated with the pathomics signature ( Figure 5A ). CRC samples with low pathomics signature scores exhibited significant enrichment in pathways related to DNA repair, proliferation, metabolism, and immune functions. Conversely, samples with high pathomics signature scores showed the activation of canonical oncogenic signaling and invasion-related pathways.

GSVA further confirmed the significant functional differences between the high- and low-pathomics signature groups ( Figures 5B, C ). Angiogenesis, epithelial-mesenchymal transition (EMT), and other invasion-associated pathways were significantly upregulated in the high-pathomics signature group. In contrast, pathways such as spermatogenesis, E2F targets, and G2M checkpoint were more active in the low-pathomics signature group.

To identify gene co-expression patterns linked to the pathomics signature, WGCNA was performed using the top 5,000 most variable genes defined by the median absolute deviation (MAD). A cluster dendrogram was generated with an optimal soft threshold power of 14 ( Supplementary Figure S11A ), resulting in 32 distinct colored modules ( Supplementary Figure S11B ). Unassigned genes were grouped into the grey module and excluded from further analysis. Correlation analysis between module eigengenes and the pathomics signature identified two modules, blue and midnight blue, as significantly associated (|bicor| > 0.1 and P-value < 0.05) ( Figure 5D ). Within these modules, gene significance was strongly correlated with module membership ( Supplementary Figure S11C ), suggesting that these genes may play pivotal roles in shaping essential biological roles related to the pathomics signature.

Subsequent functional enrichment analysis of these modules revealed distinct biological profiles ( Figure 5E ). Genes in the midnight blue module were predominantly enriched in proliferation-related pathways, whereas genes in the blue module were associated with invasion, metastasis, and immune-related processes. These findings indicate that the pathomics signature accurately reflects the underlying biological features associated with the multiple crucial hallmarks of CRC.

The MCPcounter algorithm was applied to estimate the relative abundance of stromal and immune cell subsets in relation to the pathomics signature. Both fibroblasts and endothelial cells were positively correlated with the pathomics signature and were significantly enriched in the high-pathomics signature group ( Figures 6A, B ).

To further characterize tumor microenvironment heterogeneity, the xCell algorithm was used ( Figure 6C ). Consistent with MCPcounter results, a positive correlation was observed between the pathomics signature and various endothelial and stromal cell types. Moreover, a higher proportion of lymphoid and myeloid cells was detected in samples with elevated pathomics scores. In contrast, several stem, stromal, and lymphoid cells exhibited negative associations with the pathomics signature, underscoring their relevance in reflecting CRC tumor microenvironment ( Figure 6D ).

Finally, we explored the prognostic relevance of these pathomics signature-related cells in CRC ( Supplementary Table S8 ). Survival analysis showed that cell types positively correlated with the pathomics signature were linked to poorer outcomes, whereas those negatively correlated were associated with favorable prognosis. These findings suggest that the pathomics signature captures distinct non-tumor cellular components that contribute to differential clinical trajectories in CRC.