We searched the electronic health records of two medical centers. In the first center (Hospital-1, Shanghai General Hospital – North [city center]), 100 patients admitted from July 2017 to June 2018 were randomly selected as the training cohort, and 50 patients from July 2018 to December 2018 were randomly selected as the internal validation cohort. In the second center (Hospital-2, Shanghai General Hospital – South [Songjiang new city]), 50 patients from July 2016 to December 2018 were randomly selected for external testing. The inclusion criteria were: 1) patients whose lung tumor tissue was completely resected by surgical operation; 2) lung adenocarcinoma diagnosed based on hematoxylin-eosin gross pathological specimen and immunohistochemical staining; 3) ≥ 2-year follow-up to obtain PFS results (or shorter in case of earlier progression); 4) with pre-surgical thin-slice (<1mm) contrast-enhanced CT scanning; 5) ≤ 2 months interval between CT scan and surgical resection.

The collected baseline data were age, gender, smoking status (non-smoker or smoker), 8^th edition cTNM stage, lung nodule type (solid or subsolid), histological subtype, and post-surgical treatment. The cTNM stage was determined by whole-body positron emission tomography (PET)-CT or whole-body CT for the head, neck, chest and abdomen, except the lower extremities. The histological subtype of adenocarcinoma was defined by the 2015 World Health Organization classification of lung cancer (14). The 150 patients in Hospital-1 underwent molecular testing for EGFR mutations (expression of exon-18, -19, -20, and -21) using the human gene mutation detection kit (Aide Biomedical Technology).

The study endpoint was PFS, defined as the period from tumor resection to tumor recurrence observed by PET-CT or whole-body CT, or death for any reason. The PFS results were obtained from patient medical records, including chest CT scans, clinical outcomes, or death records.

The institutional review boards approved this retrospective study and waived the requirement for patient informed consent in the two hospitals. Figure 1 shows the patient selection flow diagram. Figure 2 shows the study workflow.

Various CT scanners of different type were used in clinical practice, including three CT systems (Somatom Force, Siemens Healthineers; Revolution and HD750, GE Healthcare) in Hospital-1, and two CT systems (Somatom Flash, Siemens Healthineers; Revolution, GE Healthcare) in Hospital-2. All study patients underwent contrast-enhanced chest CT scanning after injection of 60-80 mL contrast-media (Iopamiro 300, Bracco) into the antecubital vein at 3-4 mL/sec. The image slice thickness was 0.6 mm or 0.625 mm. Supplementary Table 1 details the acquisition protocol and reconstruction parameters.

A dedicated software package (Radiomics v1.2.3, Siemens Healthineers) running on a research platform (SyngoVia VB10, Research Frontier, Siemens Healthineers) was used to perform three-dimensional tumor segmentation and feature extraction. This package was developed based on the PyRadiomics library (http://www.radiomics.io) and subjected to the Image Biomarker Standardization Initiative (15). One radiologist with 6 years of experience in thoracic imaging semi-automatically segmented the tumor lesions on CT images by finding the lesion and clicking on it, blinded to the prognostic results. Then the software automatically extracted a total of 1,672 radiomic features for each lesion, including first-order, size and shape, and texture features. Supplementary Table 2 describes these features.

Radiomic features with good stability and repeatability can be selected as candidates for further feature selection. To select these candidate features, 50 patients were randomly chosen from the training cohort. Two radiologists with 20 and 6 years of experience independently segmented the tumor lesions of these 50 patients. Spearman’s rank correlation coefficients between the 1,672 features measured by the two radiologists were calculated to filter unstable features (16). The radiomic features with a Spearman’s r > 0.8 were considered stable and reproducible, and were chosen to build the predictive model (17).

A machine-learning method, the random survival forest was used to establish the R-model and generate a Rad-score. Random survival forest provides high predictive accuracy with nonlinear regression. It is a robust algorithm with an automated feature screening process suitable for right-censored correlated complex survival datasets compared with the commonly used least-absolute shrinkage and selection operator (LASSO)-Cox regression (18). Random survival forest is suitable for integrating high-dimensional features like radiomics features for survival analysis and risk stratification. The random survival forest model was recently used to identify risk factors and generate radiomic signatures for different diseases (19). The robustness of the R-model was validated by a five-fold cross-validation approach for tuning the optimal hyperparameters. Feature importance was evaluated by the magnitude of the log-rank test statistic by permutation (20). Accordingly, the Rad-score (range from 0 to 1) was computed and generated by the random survival forest model to represent the average of expected number of events across all trees in the established forest (21). A higher Rad-score represented a higher risk of progression or shorter PFS.

A Kaplan-Meier curve with log-rank test and univariate Cox proportional-hazards model was applied to select candidate influential factors (clinical factors and EGFR mutation status) and establish a C-model. Factors with a p-value < 0.1 were included in the multivariate Cox proportional-hazards regression to establish the C-model. Prognosis prediction based on TNM staging was evaluated.

To establish the Combi-model and avoid multicollinearity and potential correlation between the Rad-score and clinical factors, a multivariate Cox proportional-hazards model with stepwise selection using minimum Akaike information criterion was implemented to integrate the Rad-score and high-relevance clinical risk factors.

Subgroup stratified analysis was implemented to verify that the radiomics model can independently predict prognosis in patients with solid and subsolid tumors. The patients in each group were classified by Rad-score as high or low risk of disease progression. The prognosis difference between risk levels was assessed in each group by Kaplan-Meier survival analysis with the log-rank test.

Clinical characteristics between datasets were compared using independent sample t-test or Wilcoxon test for continuous variables depending on the normality test, whereas Chi-square or Fisher’s exact test was used for categorical variables.

The R-score’s predictive ability was evaluated by Harrell’s concordance index (C-index) in the training, internal validation, and external test cohorts. The risk stratification capability was assessed by using the Kaplan-Meier survival curve and log-rank test. The optimal R-score cut-off values were determined by X-tile software (version 3.6.1, Yale University) (22) in the training cohort and applied to the internal validation and external test cohorts.

The proportional-hazards assumption was first verified by the Schoenfeld residuals test to evaluate model performance. Time-dependent receiver operating characteristic (ROC) curves, the integral of the area under the curves (iAUC) at various timepoints and C-index were generated and calculated in the training and internal validation cohorts to compare the prognosis discriminative ability of the three models (R-, C- and Combi-). The goodness-fit of all the three models was illustrated by the calibration curve by calculating the actual and predicted probabilities of PFS at 2 year. Decision curve analysis (DCA) was used to evaluate and compare the three models incorporating clinical net benefits. Finally, a nomogram of the Combi-model was built to visually represent the final predictive model making it convenient for clinicians to identify new patients.

Statistical analysis was performed with Python Scikit-survival v0.13.2 (https://scikit-survival.readthedocs.io), Lifelines library v0.25.5 (https://lifelines.readthedocs.io), and R package v3.6.0 (http://www.r-project.org). Detailed descriptions of the software packages and functions are listed in the Supplementary Methods. All statistical tests were two-sided, and statistical significance was set at p<0.05.

The training, internal validation and external test groups comprised 150 (61.8 ± 8.6 years old), 50 (60.8 ± 9.7), and 50 (59.9 ± 10.3) patients, respectively. All patients had histologically proven lung adenocarcinoma. Whole-body PET-CT was performed for 51/100 (51%) training, 23/50 (46%) internal validation, and 35/50 (70%) external test patients. The remaining patients underwent whole-body CT except for the lower extremities. At two years, 69/100 (69%) training, 37/50 (74%) internal validation, and 36/50 (72%) external test patients were progression-free. PFS was 409.1 ± 202.1 days, 413.8 ± 188.8 days and 541.8 ± 364.7 days, respectively. Table 1 shows the association between patient characteristics and PFS.

The feature stability analysis yielded 597 radiomic features out of 1,672 with a Spearman’s r > 0.8 ( Supplementary Table 2 ), that were considered stable features indicating good interobserver consistency. From these stable features, the 20 features highly contributive to PFS were selected to create an R-score using a random survival forest algorithm ( Supplementary Figure 3 ). The optimal R-score cut-off value to discriminate high and low progression risk was set to 3.684 based on the training cohort determined by X-tile (22), and the patients were subsequently stratified as high or low risk for disease progression ( Supplementary Figure 4 ). Kaplan-Meier analysis showed that this cut-off value could discriminate high- and low-risk patients in the training (p<0.001), internal validationt (p=0.0014), and external test (p=0.045) cohorts.

In the training cohort, the R-model reached an average AUC of 0.82 for discriminating patients with high- and low-risk for disease progression ( Figure 3A ), and an iAUC of 0.78 (95%CI: 0.67 to 0.87) and 0.92 (0.77 to 1.0) at one and two years, respectively ( Figure 3B ). In the internal validation cohort, R-model reached an average AUC of 0.77, and iAUC of 0.86 (0.74 to 0.99) and 0.70 (0.41 to 0.98) at one and two years, respectively ( Figure 3C ). In the external test cohort, the R-model showed an iAUC of 0.92 (0.70 to 0.98) and 0.90 (0.69 to 1.0) at one and two years, respectively ( Figure 3D ).

For the training and external test cohorts, the C-index was 0.781 and 0.778, respectively. All the C-indexes were higher than 0.7, indicating good predictive performance. The calibration curves showed visually good agreement between the predicted and actual labels in the training, internal validation and external test cohorts ( Supplementary Figure 5 ).

Kaplan-Meier survival analysis showed that lung nodule type and cM stage were significantly associated with PFS. The PFS was significantly shorter for patients with solid tumors than those with subsolid tumors (p=0.002), and patients with distant metastasis than those without (p=0.011). Supplementary Figure 1 shows the corresponding Kaplan-Meier curves.

The univariate Cox regression analysis showed that sex, lung nodule type (subsolid or solid), cT, cN, and cM stage were significantly associated with PFS (all p<0.1) ( Supplementary Table 3 ), so they were included for multivariate regression analysis to establish the C-model. The C-index based on TNM stage was 0.707 and 0.606 for the training and external test cohorts, respectively ( Supplementary Figure 6 ).

The global Schoenfeld test (p=0.645) showed that the C-model met the proportional hazard assumption requirement ( Supplementary Figure 7 ). Figure 4A shows the forest plots of the C-model. The calibration curves showed visually good agreement between the predicted and actual labels in the training and validation cohorts ( Supplementary Figure 8 ). The C-indexes were 0.755 and 0.739 for the training and internal validation cohorts, respectively. In the multivariate regression analysis, lung nodule type (HR=4.23, p=0.009) and cM stage (HR=4.57, p<0.001) were significantly correlated with PFS, thus were included in the Combi-model.

Schoenfeld individual tests also showed the genetic factors model adhered to the proportional hazard assumption requirement ( Supplementary Figure 9 ). The C-indexes of EGFR mutations were 0.612 and 0.528 for the training and validation cohorts, respectively. Because of these insignificant results, the genetic factor was not integrated into the Combi-model.

The Combi-model and nomogram ( Figure 5A ) were established based on Rad-score and the clinically relevant factors (cM stage, lung nodule type and histological type). A multivariate Cox regression model was used to evaluate the significance of Rad-score and these clinical factors in predicting PFS. Forest plots ( Figure 4B ) of the Combi-model illustrated that Rad-score (HR=0.034, p=0.001), cM stage (HR=3.822, p=0.002) and histological type (HR=0.191, p=0.002) were independent factors for PFS. Schoenfeld individual tests showed that the Combi-model fitted the proportional hazard assumption requirement ( Supplementary Figure 10 ). The calibration curve showed that the calibration effect of the training and validation cohorts was very good ( Supplementary Figure 11 ). The C-indexes of the Combi-model for predicting PFS in the training and external test cohorts were 0.845 and 0.837, respectively.

Finally, a nomogram combining Rad-score and clinical factors (cM stage, lung nodule type and histological type) was established ( Figure 5A ). The calibration curves showed good consistency between the predicted PFS probability and the actual result at one and two years after surgery ( Figures 5B, C ).

The pairwise comparison of C-indexes among the R-, C-, and Combi-model is shown in Supplementary Table 4 . The Combi-model performed better than the R- and C-models. In the training cohort, the C-index of the Combi-model was 0.845, higher than the R- (0.781) and C-models (0.755). In the external test cohort, the C-index for the Combi-model was 0.837, significantly higher than the R- (0.778) and C-models (0.739) (all p<0.01). DCA was used to compare and visualize the clinical net benefits of the three models ( Figure 5D ) and showed that the Combi-model gained more clinical net benefits than the other two models.

In the solid and subsolid tumor groups, Kaplan-Meier survival analysis showed that the Rad-score significantly discriminated the patients into high and low risk of progression (p=0.0045 and p<0.001, respectively) ( Supplementary Figure 12 ).