The ethics committee of our institute approved the study and waived written informed consent due to the retrospective design.

We reviewed the electronic medical records of HCC patients who received combined treatment of TACE and TKI from Sep. 2015 to Dec. 2019 at our institute (Union Hospital, Tongji Medical College, Huazhong University of Science and Technology). The inclusion criteria were as follows: 1) age, ≥ 20 years; 2) tumors confined to the liver without macro-vascular invasion or extra-hepatic metastasis; 3) tumors are measurable by the modified Response Evaluation Criteria in Solid Tumours (mRECIST); 4) Eastern Cooperative Oncology Group (ECOG) performance status of 0 or 1, Child-Pugh scores ≤7 points and adequate organ function. Those without complete medical records or high-quality CT images in electronic format were excluded.

Two radiologists reviewed pre-treatment and post-treatment CT images to evaluate short-term tumor response according to the mRECIST. Any inconsistency of assessment results was resolved by consensus. Tumor response was evaluated every 8 weeks. Overall survival (OS) was defined as the time from the date of treatment to the date of death without regarding the cause of death, and censored at the date of last follow-up for survivors. The regimen of TACE plus TKI, response evaluation, clinical data and CT data acquisition are detailed in Supplementary Material.

The Region of interest (ROI) of primary tumor, defined as enhanced area in arterial phase CT images in accordance with mRECIST, was manually delineated by two experienced radiologists (XW X and QQ R) using a 3D Slicer software (A. Fedorov et al., 2012). To be consistent with deep learning features, three consecutive slices with the maximum cross-sectional area of the tumor lesion were selected. The two observers repeated the same procedures 2 weeks later and any disagreement was resolved through consultation. The brightness, the size and of the image were standardized and the noise in the image was removed using the methods reported in literature (H. Koyuncu and Ceylan, 2018). In brief, resegmentation refers to the process whereby only pixels within a specified grey value range (−1,000, 400) are retained to exclude irrelevant organs and objects. The CT images’ appropriate window wide and center were adaptively adjusted based on the tumor region’s Hounsfield unit values. The images were then subjected to imaging normalization (the intensity of the image was scaled to 0–255) to avoid data heterogeneity bias. Histogram equalization was used to improve the brightness and contrast of the image for practitioners to analyze. CT images are mainly affected by quantum noise, arising from the variability of the electronic density of tissue voxels, statistically represented by a random Gaussian process. We used Gaussian filter to remove the noise in the image. The images with informative slices (three consecutive axial slices with maximum tumor area) corresponding to the segmented tumor region were cropped to 224 mm × 224 mm using a bounding box spanning the whole tumor area.

Six commonly used pre-trained convolutional neural networks (CNNs) (Y. Hu et al., 2021; T. N. Sainath et al., 2015), including InceptionResNetV2, InceptionV3, Resnet50, VGG16, VGG19, and Xception, were pretrained on ImageNet, which contains a large number of object categories and manually annotated training images. When performing deep learning feature extraction, we treated the pre-trained network as an arbitrary feature extractor, allowing the input image to propagate forward, stopping at the pre-specified layer, and taking the outputs of that layer as our features. After removing the last fully connected layer, we got feature maps of CT images with the maximum area of the tumor lesion, which corresponded to location invariance in the input layer. After global pooling, each feature map vector was transformed to a maximal raw value. The representational deep learning features refer to a total of 2048 (Resnet50, InceptionV3, and Xception), 1,536 (InceptionResNetV2) or 512 (VGG16, VGG19) features were converted from feature maps to numeric values.

Handcrafted radiomics features were automatically computed from the radiologist-drawn ROIs using the Pyradiomics package implemented in Python. Defined radiomics features with or without wavelet filtration were extracted in accordance with feature definitions described by the image biomarker standardization initiative (IBSI) reporting guidelines (A. Zwanenburg et al., 2020). Features were divided into three groups: (I) first-order statistics; (II) shape features; and (III) second-order features: gray level co-occurrence matrix (GLCM), gray level run length matrix (GLRLM), gray level size zone matrix (GLSZM), gray level dependence matrix (GLDM), neighborhood gray tone difference matrix (NGTDM).

The ROI images were adjusted to evaluate the impact of perturbations on feature robustness. We tested three perturbations as follows: 1) slice thickness (S): CT images were reconstructed contiguously at 1, 2, 3 and 5 mm section thicknesses; 2) rotation (R): The image and mask were rotated in the axial (x, y) plane, over a set angle θ [−30°, −15°, 15°, and 30°]; 3) segmentation (Seg): ROIs were automatically expanded or shrinked by 20% (A. Zwanenburg et al., 2019).

The Intra-class Correlation Coefficient ICC was chosen to ensure absolute agreement and not only consistency across perturbations. According to the guidelines (T. K. Koo and Li, 2016), all the features with an ICC of more than 0.85 for all tested perturbations were selected for further study analysis. Raw feature vectors were further standardized by centering on the mean and scaling to unit variance.

To further reduce feature dimension, the following steps were performed: 1) removing robust features with zero median absolute deviation (MAD); 2) only considering the top 20% features selected by univariate analysis; 3) algorithm-based feature selection; 4) the wrapper feature selection method based on the recursive feature addition algorithm to select the most predictive features. The features were fed to machine learning classifiers and the performance was evaluated by the area under the receiver operating characteristic curve (AUC). A 10-fold cross validation was used in the feature dimension step to avoid data leakage and overestimation.

The algorithm-based feature selectors included ReliefF (RELF), Fischer Score (FSCR), Gini index (GINI), Chisquare score (CHSQ), joint mutual information (JMI), conditional infomax feature extraction (CIFE), double input symmetric relevance (DISR), mutual information maximization (MIM), conditional mutual information maximization (CMIM), interaction capping (ICAP), t-test score (TSCR, only for binary classification), minimum redundancy maximum relevance (MRMR), and mutual information feature selection (MIFS). These selectors take a filter-method approach for feature selection. The filter method filters out the irrelevant feature and redundant columns from the model by using different metrics through ranking.

Twelve supervised machine learning classifiers, including Nearest Neighbors, Support Vector Classifiers (SVC) with linear or radial basis function (RBF) kernels, Gaussian processes, decision trees, random forests, multilayer perceptrons, AdaBoost, naïve Bayes, quadratic discriminant analysis (QDA), XGBoost, and logistic regression, were then used to train models for predicting short-term disease control. These classifiers were all imported from scikit-learn implemented in Python (version 3.6.4) (A. Abraham et al., 2014). During the model debugging, samples were shuffled to ensure data randomization. We adopted the Synthetic minority over-sampling technique (SMOTE) (N. V. Chawla et al., 2002), one of the commonly-used oversampling algorithms, to achieve class balance during the cross-validation step.

The terminology of each predictive model was consistent with its feature exactor, selector, and classifier. For example, VGG19_FSCR_QDA was a model trained by the QDA classifier, with features selected by FSCR and exacted by VGG19. The predictive performance of the models and their stability was evaluated by the AUC and relative standard deviation (RSD), respectively. RSD was calculated according to the formula: RSD = (sd_AUC/mean_AUC) × 100, where sd_AUC and mean_AUC were the standard deviation and mean of the ten cross-validated AUC values, respectively. Accuracy, sensitivity, specificity, precision, and F1 score were also calculated to further evaluate the selected model (Sokolova and Japkowicz, 2006).

Continuous variables with normal distribution were presented as mean ± SD (standard deviation) and those with abnormal distribution were presented as median (range). The continuous variables were compared using the t test or Kruskal-Wallis tests. Non-continuous variables were compared using the Pearson X² test or Fisher’s exact test.

Survival curves were plotted using the Kaplan-Meier method and compared using the log-rank test. Cox proportional hazard analysis was used to identify factors associated with survival. A p-value of less than 0.05 was considered statistically significant.

A total of 103 HCC patients (92 males and 11 females; age (mean ± SD): 52 ± 9 years) who received combined treatment of TACE and TKI were enrolled in this study. Of these, 72 were identified as disease control (1complete tumor response, 54 partial tumor response, and 17 stable diseases) based on mRECIST, yielding a disease control rate (DCR) of 69.9%. The rest were identified as progressed disease (PD, 30.1%). Clinical and tumor characteristics for all patients are listed in Table 1. The clinical and tumor characteristics differences between PD and non-PD groups are statistically insignificant.

Our study setup consists of three parts: 1) feature extraction and robustness analysis; 2) constructing models for predicting disease control, performance analysis, and identification of optimal models; 3) OS prediction performance analysis using the optimal models. Figure 1 shows the workflow of the study.

A consistency test was applied to evaluate feature robustness. Imaging perturbations produced a slight impact on the stability of radiomics features, with ICC of 0.93 ± 0.11 for S, 0.94 ± 0.15 for R, and 0.96 ± 0.22 for Seg, respectively. High stability was also observed in S and R perturbations for deep learning features extracted with Rnest50, with ICC of 0.89 ± 0.09 and 0.86 ± 0.12, respectively. However, Seg perturbations had moderate impact on the stability of deep learning features extracted from Rnest50, with an ICC of 0.80 ± 0.14. The results of robustness evaluation for features from all extractors were summarized in Supplementary Table S1.

Figure 2 shows the results of feature robustness analysis with ICC cutoff value of 0.85. There were 718/851 (84.37%) robust features in radiomics group. In deep learning group, the highest percentage robust features is 38.87% (199/512) from VGG19 by using the same cutoff value of ICC, followed by 35.11% (719/2048) from Resnet50, 34.38% (176/512) from VGG16, 30.62% (627/2048) from Xception, 13.61% (209/1,536) from InceptionResNetV2, and 7.57% (155/2048) from InceptionV3. The results of features robustness analysis with other ICCs were presented in Supplementary Figure S1. These results indicated that radiomics features were more stable than deep learning features; in addition, segmentation perturbation (Seg) seemed produce greater impact on stability in deep learning features.

A total of 156 radiomics features-based models and 936 deep learning features-based models were constructed, and those classified by the k Nearest Neighbors have excellent performance for predicting short-term disease control, reached a median value of AUC of 0.85 (range: 0.64–0.94) (Supplementary Figure S2) and median RSD of 1.87 (range: 0.26–11.31).

The Radiomics_GINI_Nearest Neighbors (RGNN) was identified as optimal model in radiomics group, with a cross-validated AUC of 0.87, RSD 1.30, accuracy 0.89, sensitivity 0.88, specificity 0.90, precision 0.87, and F1 score 0.89. (Figures 3A,B). Radiomics_JMI_Nearest Neighbors had a better AUC value of 0.88, but a higher RSD value of 3.59. The Resnet50_MIM_Nearest Neighbors (RMNN) was identified as the optimal model in deep learning group, with a cross-validated AUC of 0.94, RSD 0.26, accuracy 0.92, sensitivity 0.97, specificity 0.90, precision 0.83, and F1score 0.92 (Figures 3C,D). The Resnet50_JMI_Nearest Neighbors had a comparable AUC value 0.94, but a higher RSD value of 0.86.

The list of all feature selectors was in Supplementary Table S2; the ML methods’ parameter settings and tuning range were presented in Supplementary Material. The predictive performance of models constructed by other combinations of CNNs, selectors, and classifiers was Supplementary Figure S3.

For 99 patients with survival data, the median follow-up time was 15 months (range: 10–24 months). The results of the Kaplan-Meier survival analysis are presented in Figures 4A,B. There was a statistically significant survival advantage for Radiomics_GINI_Nearest Neighbors (p = 0.006) and Resnet50_MIM_ Nearest Neighbors (p = 0.033). Cox proportional hazard analysis showed that Radiomics_GINI_Nearest Neighbors (HR, 2.49; 95% CI, 1.36–4.55; p = 0.003) and Resnet50_MIM_Nearest Neighbors (HR, 1.83; 95% CI, 1.05–3.17; p = 0.032) was independently associated with overall survival (Supplementary Table S3).