Institutional review board approval was obtained for this retrospective study, with a waiver for the informed consent requirement. A total of 493 consecutive pathologically confirmed advanced stage (IIIB–IV) lung adenocarcinoma patients were recruited retrospectively from January 2014 to December 2017. The inclusion criteria were as follows: (1) age >18; (2) Eastern Cooperative Oncology Group (ECOG) performance status of 0–2; and (3) restricted therapeutic regimens: patients with TKI-sensitive epidermal growth factor receptor (EGFR) mutations or anaplastic lymphoma kinase (ALK) rearrangement accepted TKI therapy initially, and the rest of the patients accepted platinum-based chemotherapy initially. The exclusion criteria were as follows: (1) examination by unassigned CT scanners (n = 127); (2) previous anticancer therapy (n = 25); (3) incomplete clinical data (n = 92); (4) difficulty in distinguishing boundary of regions of interest (ROIs) (n = 51); and (5) loss of follow-up (n = 33). Ultimately, 165 patients were included in this study with no ALK rearrangement patients ( Figure S1 ). The clinical data collected for analysis included sex, age, ECOG, TNM stage, smoking status, TKI-sensitive EGFR mutations, tumor diameter, location, margin, lobulation, spiculation, air-bronchogram, pleural invasion, lymph node metastases, brain metastases, liver metastases, and bone metastases. The patients were randomly divided into two individual cohorts for training and validation at a ratio of 7:3 through computer-generated random numbers. The workflow of the radiomic analysis is illustrated in Figure 1 .

Contrast-enhanced CT images were acquired from Toshiba Aquilion One, Toshiba Aquilion 64 (Toshiba Medical Systems) or Phillips Brilliance iCT 256 (Philips Medical Systems) scanners. The scanning parameters were as follows: 120 kVp; 100–200 mAs; detector collimation of 64, 256, or 320 × 0.625 mm; field of view of 350 × 350 mm; matrix of 512 × 512 and reconstructed slice thickness of 2 mm. Contrast-enhanced CT scanning was performed with a 25-s delay after the injection of 85 ml of non-ionic iodinated contrast material (350 mg iodine/ml, Omnipaque, GE Healthcare). All images were exported to the Picture Archiving and Communication System (PACS) workstation (IMPAX, AGFA).

Image pre-processing was performed to enhance feature robustness and reduce feature dependence on scanner variations. Each voxel corresponded to a volume of 1.0 mm * 1.0 mm * 1.0 mm with a linear interpolation algorithm, then a Gaussian filter was used to remove noise. The gray level was consistent across the different scanners; therefore, gray level normalization was not required here.

Three-dimensional (3D) contours of the tumor regions of interest (ROIs) were delineated manually in reference to pulmonary and mediastinum windows (window width and window level of 1,500 and −450 HU on pulmonary window, while the window width and window level of 400 and 40 HU on mediastinal window). Segmentation was strictly performed by two chest radiologists (W.XT. with 7 years of experience and H.D. with 13 years of experience in chest CT) who were blinded to all patients’ information. The radiologists delineated the boundaries of the tumors on a transversal plane using itk-SNAP (version 3.4.0, www.itk-snap.org) software ( Figure S2 ). The image biomarker standardization initiative (IBSI) was regarded as reference and taken into consideration in most of the data processing, images feature, and biomarker selection procedure.

A total of 396 radiomic features were generated automatically using in-house software (Artificial Intelligence Kit, A.K., GE Healthcare) from ROIs. Features were classified into the following three categories: (a) morphological features (n = 9); (b) first-order features (intensity features, n = 42); and (c) texture features (n = 345). The details are given in Figure S3 .

Inter-/intra-class correlation coefficients (ICCs) were used to evaluate the inter-observer and intra-observer agreement. To assess inter-observer reproducibility, the ROIs of 30 randomly chosen images were performed by the two chest radiologists independently; to evaluate intra-observer reproducibility, they repeated the same procedure at an one-month interval. An ICC >0.75 was considered as good agreement. Stable and reproducible features were entered in the subsequent analysis.

Least absolute shrinkage and selection operator (LASSO) Cox regression analysis was utilized to select effective and predictable features and establish a model in the training cohort. Features with non-zero coefficients were chosen based on 10-fold cross-validation ( Figure S4 ). The radiomics signature (Rad-score), which was calculated via a linear combination of the selected features that had been weighted by their respective coefficients, represented quantitative ROI characteristics of each patient.

(1) The patients were divided into high- and low-risk subgroups in the training and validation cohorts according to the Rad-score, and the optimized threshold values were determined using X-Tile software (version 3.6.1, Yale University). Then, Kaplan–Meier OS curves and log-rank analyses were performed to assess the prognosis of subgroups. (2) The receiver operating characteristic (ROC) curves were plotted, and areas under the curves (AUCs) were calculated for predictive validity assessment of survival at 1-, 2- and 3-year time points in the training and validation cohorts. (3) The validation cohort was further divided into mutated EGFR subgroup and wild type EGFR subgroup. Harrell’s concordance index (C-index) and ROC curves for 1-, 2- and 3-year survival were used to compare the performances of radiomic model in both subgroups.

The clinical model was built by Cox proportional hazard regression to compare with the radiomic signature. Sex, age, ECOG, TNM stage, smoking status, TKI-sensitive EGFR mutations, tumor diameter, location, margin, lobulation, spiculation, air-bronchogram, pleural invasion, lymph node metastases, brain metastases, liver metastases, and bone metastases in the training cohort were first analyzed by univariate Cox regression. Only significant factors (p < 0.05) from univariate Cox regression were entered into the multivariate Cox regression analysis. AUCs were calculated for the clinical model in the same way and compared with that of radiomic signature using DeLong test.

The combined model in the form of nomogram for 1-, 2- and 3-year overall survival rate predictions was generated on the basis of the Rad-score and the clinical factors with P <0.05 in univariate Cox regression. Backward multivariate cox regression was used again, and the factors with P <0.05 were incorporated into the nomogram. The discriminative power of the predictive model was evaluated by C-index with 95% confidence intervals and AUCs in both cohorts. The calibration curves were plotted to explore the calibration degree of the combined model for the 1-, 2- and 3-year OS rates.

All statistical analyses were performed using R statistical software version 3.6.2. The “glmnet” package was used for executing the LASSO Cox algorithm. For the baseline characteristic analyses, the normality of data was assessed by the Shapiro–Wilk test. Differences between the training and validation cohorts were assessed by using independent-sample t-test and chi-square test, where appropriate. Performances of the models were evaluated by C-index. A two-sided p-value <0.05 was considered statistically significant for all comparisons.

A total of 165 patients were enrolled in the study. The cohort consisted of 94 men and 71 women with a mean age of 58.1 years (range of 34–78 years). The longest follow-up period was 72 months, and the mean was 19.2 months. The training cohort included 115 patients; the validation cohort included 50 patients. Patients’ clinical characteristics are reported in Table 1 . There were no significant differences except ECOG performance status in clinical factors between the two cohorts.

After assessing the reproducibility, 324 features with both inter- and intra-class correlation coefficients >0.75 were reserved. In the training cohort, 13 features were evaluated to construct a radiomic signature through the LASSO Cox algorithm: [1] Range; [2] skewness; [3] GLCMEntropy_ AllDirection_offset1_SD; [4] GLCMEntropy_angle135_offset1; [5] Correlation_AllDirection_offset4_SD; [6] GLCMEnergy_angle45_offset7; [7] GLCMEntropy_angle45_offset7; [8] sumAverage; [9] ShortRunLowGreyLevelEmphasis_AllDirection_offset1_SD; [10] ShortRunEmphasis-_AllDirection_offset4_SD; [11] ShortRunHighGreyLevelEmphasis_AllDirection_offset7_SD; [12] ShortRunHighGreyLevelEmphasis_angle0_offset7; and [13] Sphericity. The formula of Rad-score is illustrated in Table S1 .

The C-indexes of the radiomic signature were 0.746 (95% CI, 0.699 to 0.792) in the training cohort and 0.677 (95% CI, 0.597 to 0.766) in the validation cohort.

(1) The patients were classified into high- and low-risk subgroups according to the Rad-score at a cut-off point of 0.15 according to X-Tile. Kaplan–Meier analysis revealed significantly different subgroup OS in both training cohort (p < 0.0001, log-rank test) and validation cohort (p < 0.0001, log-rank test), as shown in Figure 2 .

(2) The ROC curves of the two cohorts for 1-, 2- and 3-year survival are plotted in Figure 3A .

(3) The C-index for mutated EGFR subgroup was 0.629 (95% CI, 0.476 to 0.782) in validation cohort; the C-index for wild type EGFR subgroup was 0.662 (95% CI, 0.537 to 0.787) in validation cohort. The AUCs of the subgroups for 1-, 2- and 3-year survival are shown in Table S2 . Although there was no significant difference by the DeLong test in all pairs, the AUCs of wild type EGFR subgroup were superior to that of mutated EGFR subgroup for all cut-off time points.

The variables with p values <0.05 in the univariate analysis, namely, ECOG, TKI-sensitive EGFR mutations (treatment methods), pleural invasion, and brain metastases, were entered into the multivariate analysis. The multivariate Cox proportional hazard model showed that all the entered factors were identified as independent predictors of OS ( Table 2 ).

For the clinical model, the C-indexes were 0.718 (95% CI, 0.669 to 0.766) in the training cohort and 0.698 (95% CI, 0.603 to 0.792) in the validation cohort. The ROC curves for the clinical model are depicted in Figure 3B . The comparison of AUCs for radiomic signatures and the clinical model are shown in Figure 4 . In the validation cohort, although there was no significant difference by the DeLong test in all pairs, the AUC of the radiomic signature for 1-year prediction was inferior to that of the clinical model, but the 2-year prediction was superior to that of the clinical model. The prediction efficiencies of both models for 3-year survival were not satisfactory.

Clinical parameters with p <0.05 in the univariate Cox regression (ECOG, treatment methods pleural invasion and brain metastases) and Rad-score were included in the construction of the combined model using backward multivariate Cox regression ( Table 3 ). The nomogram was showed in Figure 5 .

The C-indexes of the combined model were 0.799 (95% CI, 0.757 to 0.84) in the training cohort and 0.733 (95% CI, 0.656 to 0.81) in the validation cohort, which were higher than those of the other two models. ROC curves are delineated in Figure 3C . Figure 4C shows that the combined model improved the accuracy for 2- and 3-year survival predictions compared with the radiomics and clinical models. The calibration curves of the nomogram demonstrated good consistency between predicted and observed results ( Figure 6 ).