With approval from the local ethics committee of our hospital, this retrospective study waived the requirement for informed consent. The study design is illustrated in Figure 1.

From January 2010 to December 2015, consecutive patients with histologically proven stage II CRC and available follow-up information were included. All patients received surgucal resection within 2 weeks after preoperative CT scan. The exclusion criteria were as follows: (1) patients previously treated with any anticancer therapy (n = 44), (2) patients with co-malignancy (n = 27), (3) patients with lack of available baseline demographics and CT images (n = 132), and (4) patients lost to follow-up within 3 years (n = 65) and those with incomplete clinicopathological data (n = 32).

Consequently, 299 patients (median [interquartile range {IQR}] age, 59 [30–82] years) were enrolled in our study, including 187 male (median [IQR] age, 59 [31–82] years) and 112 female (median [IQR] age, 58 [30–80]; Table 1). For temporally independent validation, patients were randomly separated into two cohorts at a 7 to 3 ratio: 210 for training and 89 for validation.

In this study, the clinicopathological risk factors included gender, age, T stage, CT-reported tumor location (Ascending colon, Transverse colon, Descending colon, Sigmoid colon, and Rectum), histologic grade of the tumor, smoking history, hypertension history, diabetes history, family history of cancer, internal obstruction, or perforation (IOP) status, number of lymph nodes examined (≥12 vs. <12), LVI status, PNI status, mismatch repair status, Ki-67 expression level, and history of postoperative adjunctive chemotherapy (present/absent). Laboratory analysis included tests for carcinoembryonic antigen (CEA) with the values of 5 ng/mL, carbohydrate antigen 724 (CA724) with 6 U/mL, carbohydrate antigen 242 (CA242) with 20 U/mL and carbohydrate antigen 199 (CA199) with 39 U/mL. The histologic grade was on the basis of the World Health Organization (WHO) classification of the digestive system tumors, 4th edition (21). Ki-67 was identified as positive expression when staining was equal to or above 40% of the specimen, while <40% denoted negativity (22). DNA mismatch repair (MMR) status was assessed by immunohistochemical staining for MMR gene protein products (MLH1, MSH2, MSH6, and PMS2) expression (23). The loss of an MMR protein in tumors was defined as high-frequency microsatellite instability (MSI-H), whereas intact MMR proteins in tumors was defined as low-frequency MSI (MSI-L). These datasets were obtained from the institutional archives.

After surgical resection, all patients were followed-up for at least 3 years. Patients were followed-up by recording their CEA levels and evaluating contrast CT images in the first month after surgery and every 3–6 months thereafter. The end-point was time to recurrence, which was defined as the prognostic performance of the imaging features for distant metastasis, local recurrence or an atypical finding with histopathological confirmation.

DFS was calculated from the date of surgery until either the date of confirmed clinical recurrence or time of last available contact.

The pretreatment abdominal with/without pelvic contrast-enhanced CT scans were obtained using a varied set of CT scanners (Supplementary Material 1).

For tumor segmentation, all portal venous-phase CT images with DICOM format were retrieved from the picture archiving and communication system (PACS) at our institution, because of well-differentiation between tumor tissue and adjacent normal bowel wall. Volume of interests (VOIs) were semiautomatically delineated using the open-source software 3D Slicer (version 4.9, www.slicer.org). This was performed by a board-certified radiologist (reader 1, C.X.N.) with 6 years of experience in abdominal radiology and then reviewed and modified by another experienced radiologist (reader 2, F.S.X., with 10 years of experience in abdominal diagnosis); both were blinded to the clinicopathologic and outcome details. To assess feature reproducibility, segmentation was repeated in 20 randomly selected patients by another radiologist (reader 3, S.Q., with 10 years of experience in abdominal diagnosis). On the basis of the feature extraction by reader 1 and reader 3, the interobserver ICCs were calculated. Additionally, reader 1 repeated the assessment of the same 20 randomly selected patients 1 month later to evaluate intraobserver reproducibility.

After tumor segmentation, 1561 radiomics features (RFs) per case were extracted from the region of interest (ROI) with the open-source python package “Pyradiomics V1.3.0” (http://www.radiomics.io/pyradiomics.html). The following radiomic features were included for feature computation: (1) Shape-based features; (2) First-order statistics features; (3) gray level co-occurrence matrix-based features (GLCM); (4) gray level run-length matrix-based features (GLRLM); (5) gray level size zone matrix (GLSZM); (6) neighboring gray tone difference matrix (NGTDM), and (7) gray level dependence matrix (GLDM). The detailed feature extraction methodology was shown in Supplementary Material 2.

A two-step procedure of feature selection was performed for dimensionality reduction. Firstly, we calculated the intra- and interclass correlation coefficients (ICCs) for all 1561 radiomics features based on the abovementioned resegmentations to remove the unstable features. Only features with both intra- and interobserver ICC values >0.90 were initially selected. Secondly, considering that the extracted features were high dimensional, three steps were adopted to reduce the dimension of the features: (1) Kruskal-Wallis test was first used to screen the image features with statistically significant differences (P < 0.05). The objective is to remove the irrelevant or poorly correlated characteristic parameters. (2) In consideration of the possibility of repeated expression of lesion information among features, Spearman analysis was adopted as the correlation analysis to remove image features with correlations <0.8. (3) The Least Absolute Shrinkage and Selection Operator (LASSO) logistic regression algorithm (24, 25), with penalty parameter tuning conducted by 10-fold cross validation, was used to preidentify the top-ranked features. The candidate predictive features with a zero-fit weight were selected. Thereafter, a Rad-score for each outcome was built via a linear combination of selected features and coefficient vectors. The Rad-score cutoff values for the classification of high- and low-risk groups were chosen according to the Youden index criterion (26).

All radiomics feature selection and model construction processes were performed in the training cohort and then evaluated in the validation cohort.

Chi-square test, Student's t-test, Fisher's exact test and Mann-Whitney U-test were used for categorical and continuous variables, where appropriate. For survival analyses, we used the Kaplan-Meier method to analyze DFS and the log-rank test to compare survival curves. A multivariate logistic regression analysis with Cox regression coefficients was performed to construct nomograms combining clinicopathological factors and RFs. To evaluate the predictive accuracy of the nomograms, a calibration curve with bootstrapping, which measured the agreement between the nomogram-predicted outcome probability and the average actual probability, was plotted. The performance of the clinical model, radiomics signature and combined model were evaluated by diagnostic accuracy, sensitivity, specificity, and the time-dependent receiver operating characteristic (ROC) curve. Differences between various AUCs were compared with the DeLong test. All statistical analyses were performed using R software (version 6.1, www.r-project.org). A two-tailed P < 0.05 was considered a statistically significant difference. Due to sample imbalance (the number of recurrences was much smaller than that of nonrecurrences), the synthetic minority over-sampling technique (SMOTE) (27, 28) was used to overcome classification bias in favor of the majority class. Further details are included in Supplementary Material 3.

The detailed clinicopathological characteristics and demographics of patients in the training and validation cohorts are shown in Table 1. The median DFS time was 1,024 days (range, 23–1,978 days), and 81 (27.1%) of 299 patients experienced recurrence during the follow-up period.

After the feature robustness analysis, we adopted the 118 most stable features with both intra- and interobserver ICC values >0.90 for further analysis. Then, the 114 most significant features selected by the LASSO logistic regression analysis were adopted in the Rad-score calculation formula (Figure 2, Supplementary Table 1). The Rad-score calculation formulas are listed in Supplementary Material 4.

Accordingly, Rad-score of 0.374 was calculated to differentiate between high-risk and low-risk group. Ninety-one (43.3%) patients were assigned to the high-risk group, and among them, 61 (67.0%) had developed recurrence by the overall endpoint. A total of 119 (56.7%) patients were assigned to the low-risk group, and among them, 1 (0.8%) had developed recurrence by the overall endpoint. The recurrence risk in the high-risk group was significantly higher than the low-risk group (odds ratio 239.933, P < 0.001). We performed the same analyses in the validation cohort. The 3-year DFS was 55.9% for the high-risk group and 0% for the low-risk group (HR 130, 95% CI 18–940; P < 0·0001; Figure 3).

Univariate Cox regression analysis in the training cohort showed that the Rad-score, CA242, CA724, CA199, and CEA levels, mismatch repair status, T stage, LVI, PNI, IOP status, number of lymph nodes examined, and adjuvant chemotherapy were significant prognostic factors for DFS (Table 2). Subsequently, all the relevant factors were entered into the multivariate Cox analysis, and the Rad-score, CA724 level, mismatch repair status, and PNI were considered independent risk factors for model building (Table 2). Detailed results of the validation cohort analyses are listed in Supplementary Table 2.

It should be emphasized that after a multivariable analysis was performed to adjust for the clinicopathological risk factors, the radiomics-based classifier remained a powerful and independent factor in both the training and validation cohorts (HR 38.08; 95% CI: 14.45–100.33, P < 0.001; HR 60.72; 95% CI: 6.62–557.15, P < 0.001). When stratified by different clinicopathological risk factors, patients in the low-risk group significantly showed better 3-year DFS than those in the high-risk group (Figure 4).

Surprisingly, the multivariable cox regression analysis showed that adjuvant chemotherapy was not an independent risk predictor in the training (HR 0.497, 95% CI 0.229–1.077; P = 0.076; Table 2) and validation sets (HR 2.560, 95% CI 0.429–15.289, P = 0.303; Supplementary Table 2). However, our radiomics-based classifier indicated that patients in the high-risk group derived a greater survival benefit from adjuvant chemotherapy (HR 0.303, 95% CI 0.181–0.506, P < 0.001; Supplementary Figure 1).

Table 3 presents the performance results obtained in the training and validation cohorts for the clinical, radiomics, and radiomics plus clinical models.

Overall, in the stratification analysis, the Rad-score showed a significant discrimination of high-risk and low-risk patients with CRC in the training cohort (AUC 0.886, 95% CI 0.840–0.931) and the validation cohort (AUC 0.874, 95% CI 0.802–0.945). The combined model yielded an AUC of 0.954 (95% CI 0.930–0.978) in the training cohort and 0.906 (95% CI 0.844–0.968) in the validation cohort, and these values were significantly greater than those obtained the model of clinical parameters alone (AUC 0.756, 95% CI 0.694–0.817, P < 0.001; AUC 0.704, 95% CI 0.586–0.823, P < 0.001) (Figure 5). Comparison of the performance demonstrated no significance difference between the training and validation cohorts (P = 0.160). Thus, the radiomics classifier could add prognostic value when combined with clinicopathological prognostic features.

To provide the physicians with a simple and quantitative approach to predict the disease recurrence probabilities for each individual patient, we developed a combined nomogram that integrated both the Rad-score and clinicopathological risk factors (Figure 6A). Notably, for the Rad-score, the variable with the largest coefficient absolute value was set as a reference, whose scale ranged was from 0 to 100. The calibration curve of the combined nomogram showed good agreement between the nomogram prediction and actual observation (Figures 6B,C). The nomogram was able to accurately predict the 3-year DFS, with a C-index of 0.872.