Shanxi Medical University Review Board approved this retrospective study. The entire cohort was obtained from the Institutional Picture Archiving and Communication System (PACS) at Shanxi Cancer Hospital from February 2016 to October 2018. The patients who had histologically confirmed ESCC (TNM stage: I-III) and underwent surgery after diagnosis, underwent pretreatment CT scans from neck to abdomen and signed informed consent. All methods were carried out following the relevant guidelines and regulations.

To determine the patients that could be included, we developed the following criteria: (1) pathologically confirmed ESCC; (2) underwent surgery for ESCC; (3) standard contrast-enhanced CT was performed preoperatively; and (4) complete clinical and follow-up information was available. We randomly divided the patients into training and validation cohorts by a ratio of about 3:1. We trained models in the training cohort and validated them in the validation cohort.

Clinical characteristics including age, gender, tumor location (upper, middle, lower), drinking history, smoking history, genetic alterations, and pathologic features including depth of invasion, TNM stage, and lymph node metastasis information were collected from patient records. These clinicopathologic characteristics are presented in Table 1.

All patients were followed up every 1–3 months during the first 2 years, every 6 months in years 2–5, and annually after that. To provide an efficient tool, which would allow earlier personalized treatment, we chose PFS as the endpoint (Sargent et al., 2005). We defined PFS from the first day of treatment to the date of disease progression (locoregional recurrences or distant metastases), death from any cause, or the date of the last follow-up visit (censored). The minimum follow-up time to ascertain the PFS was 2 months.

All patients were underwent the contrast-enhanced CT using a 64-channel multi-detector CT scanner (LightSpeed VCT, GE Medical Systems, Milwaukee, Wis, USA). The acquisition parameters were: 120 kV; 160 mA; 0.5-s rotation time; detector collimation: 64 × 0.625 mm; field of view: 350 × 350 mm; and matrix: 512 × 512. After routine non-enhanced CT, contrast-enhanced CT was performed after a 25-s delay following intravenous administration of 85 mL of iodinated contrast material (Ultravist 370; Bayer Schering Pharma, Berlin, Germany) at a rate of 3.0 mL/s with a pump injector (Ulrich CT Plus 150, Ulrich Medical, Ulm, Germany). All images were reconstructed with a thick slice of 5.0 mm. We converted the image format from DICOM to NII for feature selection without any preprocessing.

Note that segmentation is required before the extraction of quantitative radiomics features; we performed three-dimensional manual segmentation using 3D-Slicer software (https://www.slicer.org/), an open platform for medical image processing. The chief physician of Shanxi Cancer Hospital, with more than 5 years' experience in interpreting chest radiology, outlined the tumor regions for each CT image layer. The tumor segmentation was guided and verified by the specialist. The region of interest (ROI) covered the whole tumor mass, was delineated on each CT slice, and used in subsequent feature extraction.

We performed the calculation through our homemade Python scripts (Python3.6, https://www.python.org) for radiomics feature extraction based on the segmentation results. A total of 954 features were obtained by calling feature calculation in pyradiomics package (open-source python package; https://pyradiomics.readthedocs.io/en/latest/), which included the following fourcategories: (1) first-order statistics features; (2) size- and shape-based features; (3) texture features; and 4) wavelet features; and five typical matrixes: Gray-Level Co-occurrence Matrix (GLCM), Gray Level Run Length Matrix (GLRLM), Gray Level Size Zone Matrix (GLSZM), Gray Level Dependence Matrix (GLDM) and Neighboring Gray Tone Difference Matrix (NGTDM).

We built the radiomics signature with selected features in the training cohort. To reduce over-fitting or any types of bias, we applied the following two steps: First, the best features based on univariate statistical tests (2-sample t-test) between death and censoring groups in the primary cohort were selected and executed by using Matlab 2016b. Second, we used our homemade R scripts to select features that were most significant by using the least absolute shrinkage and selection operator (LASSO) method, which would be a suitable methodology for the feature selection through regression of high-dimensional data. Additionally, the accuracy of the prediction model could be improved by regularizing the features through penalized estimation. We added the L1 penalty term to the normal linear model, and the parameter lambda controls the complexity of regression. When the λ is large, it indicates no effect on the estimated regression parameters, while the λ gets smaller, most covariate coefficients were shrunk to zero. Then the remaining variables with non-zero coefficients were selected by the λ that the 10-fold cross-validation error was the most minor (Kumamaru et al., 2016; Vasquez et al., 2016). When performing 10-fold cross-validation, the training cohort was divided into 10 equal parts; each called a fold. Next, it will train a series of models. The first model was trained using the first fold as the test set and the other folds (2–10) as the training set. Then another model was constructed using 2nd fold as the test set and the1st, 3th-10th folds as the training set. The process is repeated with 3th-10th folds as a test set.

Finally, the radiomics signature was built by combining those variables in the primary cohort and validated in the validation cohort. The radiomics signature is a linear combination of selected features with respective weights, which would be calculated as a factor (Radiomics score, Rad-score) for the further prediction model. The Rad-score calculated by using the following formula:

Where F_i is the selected radiomics feature, c_i is the LASSO coefficient of F_i. Then, the assessment method of the logistic regression model is the receiver operating characteristic (ROC) curve and its area under the curve (AUC).

We calculated Rad-score for each patient with ESCC and grouped them according to two rules. (1) The patients were divided into high-risk and low-risk groups based on the median Rad-score. (2) Patients with median scores were placed in high-risk groups. The radiomics signature discriminative performance of the survival status was assessed according to the overall distribution of ESCC patients. And then, the potential association of radiomics signature and clinical feature with PFS was assessed in the training cohort and validated in the validation cohort. Kaplan–Meier survival analysis was used in these two cohorts. Stratified analyses were implemented to determine the PFS in high-risk and low-risk patient subgroups. Univariate Cox Proportional Hazards Models were performed to explore the C-index of the radiomics signature.

The nomogram with the predicting model was based on the multivariable logistic regression analysis. The following candidate factors: TNM stage (dummy variable: “0” for I, “1” for II, “2” for III), the status of clinical features, and Rad-scores were involved in a diagnostic model for preoperative prediction of ESCC. The nomogram is a graphical representation of this prediction model in the training cohort. The prognostic performance of TNM staging and clinical nomograms in the training cohort before and after the addition of the Rad-score was quantitatively measured using Harrell's concordance index (C-Index), which is commonly used for the evaluation of the discriminative power of prognostic models (Harrell, 2018). The value of the C-index could range from 0.5, which indicated no discriminative ability, to 1.0, which showed perfect ability to distinguish between the patients who suffered disease progression or death and those who did not. Bootstrap analyses with 1,000 resample were applied to obtain a C-index with a 95% confidence interval (CI) (Canty and Ripley, 2021) corrected for potential overfitting. The calibration curves were drawn for assessing the agreement between the predicted probability of 3-year PFS and actual 3-year PFS (Pencina et al., 2011).

The prognostic performance of TNM staging and clinical nomograms in the validation cohort before and after the addition of the Rad-score was tested by the above method. The calibration curve and C-index were calculated through multivariable Cox proportional hazard regression analyses. The decision curve analysis (DCA) was introduced to evaluate the quantified net benefit of our prediction model in the validation cohort (Vickers et al., 2008; Shen et al., 2018).

A heat map analysis was used to evaluate the associations between clinical data and radiomics features.

All the statistical analyses were performed using IBM SPSS software (version 26; IBM Corp, Armonk, NY, USA), Matlab 2016b (Mathworks, Natick, USA), and R software (version 4.1.1, Boston, MA, USA). In this study, 2-sample t-test was applied to confirm whether differences between death and censoring groups in the primary cohort. The differences in gender, age, TNM stage, smoking status, drinking status, location, genetic history, invasion degree and metastasis between the training and validation data sets were assessed by using the χ²-test. The following R packages were used: the glmnet package was used for the LASSO logistic regression model analysis, the pROC package was used for the ROC curves, the Hmisc package was used for the comparisons between the C-indices, the survival package was used for Kaplan–Meier survival analyses, the rms package was used for the nomograms and calibration curves, the rmda package was used to implement the DCA, and the gplots and pheatmap packages were used for heat maps.

A total of 272 consecutive patients who met the criteria (192 men and 80 women; mean age, 60.27 years ± 7.43) were included and divided into two cohorts by a ratio of 3:1 using computer-generated random numbers. A total of 204 patients were enrolled in the training cohort (146 men and 58 women; mean age, 60.22 years ± 7.28), while 68 patients were enrolled in the independent validation cohort (46 men and 22 women; mean age, 60.44 years ± 7.90). The clinical characteristics with statistics of the training and validation cohorts are summarized in Table 1. No significant differences were found between these two cohorts in terms of gender, age, history of smoking and drinking, location, genetic history, invasion degree, lymph node metastasis, and overall TNM Stage (p = 0.255–0.887). The median PFS was 35.35 months (range, 2–75 months).

A total of 954 features were extracted from CT images and might contain many redundant and highly correlated features. To find out robust and valuable features, we performed the following steps: First, 221 features were selected by univariate statistical tests (p < 0.05) (Table 2). Then, based on the LASSO logistic regression algorithm approach in the training cohort, we selected the features with non-zero coefficients. As a result, 21 radiomics features were screened out from 221 features (Table 3). The procedures of parameter tuning and feature space reduction of the regression model are illustrated in Figure 1. Then the 21 features were selected to build the radiomics signature and involved in the Rad-score-based prognostic model. The discriminative ability of the survival status based on radiomics signatures was assessed by ROC in both cohorts, respectively (Figure 2A).

Rad-score for each patient in the training cohort and the validation cohort correspondingly showed that the higher the Rad-score, the greater the probability of death (Figures 2B,C). Besides, in the training cohort, the radiomics signature from CT images yielded the highest C-index, which was 0.758 (95% CI: 0.708–0.808). In the validation cohort, the radiomics signature from CT images yielded a C-index of 0.748 (95% CI: 0.656–0.840). It showed a significant discrimination between the PFS of high-risk and low-risk patients (Figure 3).

We developed a radiomics nomogram that integrated the radiomics signature from the CT images with the traditional TNM staging system, which yielded a C-index of 0.603 (95% CI: 0.549–0.657). This nomogram significantly improved the discrimination ability in evaluating PFS (C-index: 0.770; 95% CI: 0.721–0.819) than TNM staging system (p < 0.05; Figure 4A), and showed good calibration as well (Figure 4B). Moreover, a radiomics nomogram was created by integrating the radiomics signature from the CT images with all clinical data, whose nomogram yielded a C-index of 0.680 (95% CI: 0.626–0.734). We found that the radiomics nomogram possessed good calibration and seemed to be more accurate than the clinical nomogram for evaluating PFS (C-index: 0.792; 95% CI: 0.748–0.836) with a p < 0.05 (Figures 4C,D).

In the validation cohort, the C-index of the traditional TNM staging system is 0.572 (95% CI: 0.478–0.666). We integrated the radiomics signature with the TNM staging system to produce a radiomics nomogram, which showed an improvement over the TNM staging system alone (C-index: 0.760; 95% CI: 0.673–0.847). The calibration curve of probability in PFS evaluation showed good agreement between nomogram-evaluated and actual observation (Figure not shown). While the clinical nomogram yielded a C-index of 0.605 (95% CI: 0.501–0.709) in the validation cohort and was advanced by combining with radiomics signature (C-index: 0.779; 95% CI: 0.697–0.861). The calibration curves of this nomogram showed good agreement between nomogram-evaluated and actual survival (Figure not shown). The DCA for the prediction model derived from the addition of Rad-score before and after is presented in Figure 5A. It showed that the predictive model collaborated with Rad-score had a better net benefit than that with only traditional TNM staging combined with clinical features.

The ESCC patients with similar patterns of radiomics expression were clustered through unsupervised clustering (Figure 5B). Then we organized a heat map to determine the association between radiomics features and clinical data (Figure 5B). The results showed significant correlations between signature features LHL_glcm_Idn, LHL_glszm_Large Area High Gray Level Emphasis with drinking (p < 0.001) as well as gender (p < 0.001). LHL_glcm_ldn was associated with smoking (p < 0.01). LHL_gldm_Small Dependence Low Gray Level Emphasis was associated with location (p < 0.01). LHL_glcm_Idn and LHL_glszm_Large Area High Gray Level Emphasis were associated with invasion degree (p < 0.05). HLH_glszm_Zone Entropy, LHL_glcm_Idn and LLH_firstorder_Energy were associated with overall stage (p < 0.05). In contrast, no radiomics feature was significantly associated with age, genetic history and metastasis (for all, p > 0.05). These p-values for all the correlations of radiomics features with clinical data are shown in Supplementary Table 1.