A multicenter retrospective study involving three hospitals was performed, which was approved by the ethics committee. Informed consent was waived given the retrospective nature of the study. Clinical data and ultrasound images were collected from 2,137 patients pathologically confirmed as HCC or ICC at three hospitals from January 2019 to March 2021. Among them, 226 patients (HCC = 176, ICC = 50) were included in the final analysis. The inclusion criteria were as follows: (1) being pathologically confirmed as HCC or ICC; (2) having received liver ultrasound within 1 month before surgery and the ultrasound images information being intact; (3) having not received antitumor treatments before, including liver transplantation (LT), microwave ablation (MWA), radiofrequency ablation (RFA), and transcatheter arterial chemoembolization (TACE); (4) the ultrasound images satisfying the analytical requirements and the target lesions being totally visible on the ultrasound images; and (5) no history of concurrent malignancies. The flow chart of subject inclusion and exclusion is shown in Figure 1 .

Patients from two hospitals (Henan Provincial Peoples Hospital and the First Affiliated Hospital of Zhengzhou University) (n = 187, HCC = 146, ICC = 41) were mixed and divided into training set (n = 149, HCC = 118, ICC = 31) and test set (n = 38, HCC = 28, ICC = 10) by stratified sampling at a ratio of 8:2. The patients from a third hospital (Henan Cancer Hospital) (n = 39, HCC = 30, ICC = 9) were used as an independent validation set.

The clinical data were acquired from the electronic health records, including: demographics (gender, age, history of hepatitis), laboratory tests (AFP, ALT, AST, TB, CB, and UCB) and ultrasound features (size of lesion). The laboratory tests and ultrasound imaging were examined within 1 month before surgery. The patients’ pathology information (the pathological diagnosis of HCC or ICC) was obtained from the pathology information system.

Ultrasound images of liver tumors were collected using the Color Doppler Ultrasound System with a convex array transducer (frequency range 2.5–6 MHz), including GE Logiq E9, GE Vivid E9, HI VISION Ascendus, HI ALOK ProSound A5, Philips EPIQ 5, and Aloka EZU-MT28-S1. All ultrasound scans were performed by ultrasound physicians who had over 5 years of experience in liver ultrasound. At least one original ultrasound image showing the lesion and the same image containing the measurement parameters were stored in the DICOM format.

The open-source software ITK-SNAP v.3.6.0 was used to manually delineate the region of interest (ROI) (34). First, an ultrasound physician with over 9 years of experience loaded the images into the ITK-SNAP software and manually annotated the entire lesion. Another ultrasound physician with 30 years of experience then delineated ROI in the lesions for all ultrasound images. The reproducibility of feature extraction from ROI was evaluated according to the delineation results. Both ultrasound physicians had 4 years of working experience concerning ITK-SNAP software. They were blinded to clinical history and pathology results but were aware of the purpose and design of the study. The ROI segmentation results for the representative liver lesions are shown in Figure 2 .

A researcher with 5 years of experience performed image preprocessing to eliminate variability of the ultrasound images arising from the use of different ultrasound equipment at different hospitals and to improve the reproducibility of feature extraction. First, the ultrasound images were normalized based on the mean and standard deviation. Second, the images were resampled by B-spline interpolation to 1 mm × 1 mm pixel. Finally, gray-level discretization was performed for the histogram with the bin width fixed at 25 (35).

The open-source Python package Pyradiomics v.2.1.2 was used to extract ultrasound features from each patient. The extracted features were divided into the following seven categories (36): (I) first-order statistical features; (II) two-dimensional shape features; (III) gray-level cooccurrence matrix (GLCM); (IV) gray-level run length matrix (GLRLM); (V) gray-level size-zone matrix (GLSZM); (VI) gray-level dependence matrix (GLDM); (VII) neighborhood gray-tone-difference matrix (NGTDM). Fourteen filters were applied on the original images to obtain the derivative images for each patient. Except for the shape features, the features of all categories were directly obtained from the original and derivative images. Detailed information about the feature extraction method, filters, and features is available in the Supplementary Materials 1, 2 .

Since the unit and value range varied for different extracted features, the feature values were of varying scales. To cope with this problem, we performed Z-score normalization before feature selection to ensure a relatively uniform distribution of the image features. However, all of the extracted features were high dimensional. The use of high-dimensional features might have the problems of low computational efficiency and overfitting (37). First, the features with zero variance were excluded by using the variance filtering method. Next, lasso method was performed for further dimensionality reduction of the features and the most valuable features were selected. The 10-fold cross-validation process was repeated 1,000,000 times to obtain the optimal value of parameter λ, which was introduced into the lasso method to calculate the regression coefficients of each feature. Finally, the features with nonzero coefficient were selected. The study workflow is shown in Figure 3 .

We invoked the Python scikit-learn 0.23.2 package for SVM model training and performance evaluation. The patients from two hospitals were randomly divided into training set and test set by stratified sampling at a ratio of 8:2. The patients from a third hospital were used as an independent validation set. The learning curve and the grid search were used concomitantly to select the optimized parameter combination consisting of the kernel function, coefficient of kernel function, penalty coefficient and class_weight. The specific process of parameter tuning is available in the Supplementary Material 3 .

Three models were constructed in this paper. First, the clinical model was constructed using the patients’ clinical data, including gender, age, history of hepatitis, AFP, ALT, AST, TB, CB, UCB, and the size of lesion. Second, an ultrasomics model was constructed using the ultrasomics signatures extracted and selected from the ROI delineated on the ultrasound images of the HCC or ICC patients. Finally, the combined model was built by integrating the clinical features and the ultrasomics signatures. The details of the model construction process can be found in the Supplementary Material 4 .

The three models built upon the training set were evaluated using the test set and the independent validation set. The predictive performance of the three models was evaluated by plotting the ROC and estimating the performance indicators, including AUC (along with 95% CI), accuracy, sensitivity, and specificity. An overview of the entire process is shown in Figure 3 .

SPSS 25.0 software was used for statistical analysis. The normality of continuous variables was tested using the Kolmogorov-Smirnov test. Continuous variables obeying a normal distribution were analyzed by the independent-samples t-test. Otherwise, they were analyzed by Wilcoxon’s rank-sum test. The relationships between the categorical variables were tested by using the chi-square test. The continuous variables obeying a normal distribution were expressed as mean ± standard deviation. Otherwise, the continuous variables were expressed by medians [interquartile range (IQR)]. Categorical variables were expressed as n (%). p < 0.05 indicated a significant difference. The Delong test was employed for a quantitative comparison of the ROC among the three models (38).

The reproducibility of feature extraction was evaluated using the intraclass correlation coefficient, which greater than 0.8 indicated high consistency, 0.5 to 0.79 moderate consistency, and less than 0.5 low consistency (39).

The clinicopathological features in the training set, test set, and independent validation set are shown in Table 1 . The percentages of ICC patients in the training set, test set, and independent validation set were 20.8% (31/149), 26.3% (10/38), and 23.1% (9/39), respectively. The percentages of ICC patients with a history of hepatitis were 63.8% (95/149), 73.7% (28/38), and 74.4% (29/39), respectively. The average age of patients was 57.2 ± 11.1, 58.7 ± 9.3, and 59.1 ± 11.2 in the training set, test set, and independent validation set, respectively. The three sets did not differ significantly in demographics, laboratory test results, and ultrasound features (p > 0.05).

From each patient, 1,409 features were extracted from the ultrasound images. Among the extracted features, fourteen 2D shape features were obtained only from the original images. Except for that, the features of all other six categories were obtained from one original image plus 14 derivative images. There were 18 first-order statistical features, 24 GLCM features, 16 GLRLM features, 16 GLSZM features, 14 GLDM features, and 5 NGTDM features. More information about the extracted features is listed in the Supplementary Material 2 .

Among these extracted features, 330 features with an intraclass correlation coefficient below 0.8 were first excluded. Then 16 features with zero variance were excluded using the variance filtering method. Lasso was used to reduce the dimensionality, which finally resulted in 14 features. The process of lasso feature selection is illustrated in Figure 4 , with detailed information shown in the Supplementary Figures S1–S3 .

The ROC curves of the clinical model and the ultrasomics model on the training set, test set, and independent validation set are shown in Figures 5A, B . On the test set, AUC (along with 95% CI), sensitivity, specificity, and accuracy of the clinical model and the ultrasomics model were 0.711 (0.541–0.846), 0.700, 0.714, and 0.711 vs. 0.843 (0.688–0.940), 0.900, 0.750, and 0.790, respectively. On the independent validation set, these performance indicators were 0.800 (0.641–0.911), 0.889, 0.667, and 0.718 vs. 0.730 (0.564–0.859), 0.667, 0.700, and 0.692, respectively. According to the results above, the ultrasomics-based model outperformed the clinical model in AUC and accuracy on the test set. The situation was just the opposite on the external validation set. The performance of the clinical model and the ultrasomics model in differentiation between HCC and ICC on the training set, test set, and independent validation set is shown in Table 2 .

Figure 5C displays the ROC curves of the combined model on the training set, test set, and independent validation set. The AUC (along with 95% CI), sensitivity, specificity, and accuracy of the combined model on the test set and the independent validation set were 0.936 (0.806–0.989), 0.900, 0.857, and 0.868 vs. 0.874 (0.733–0.961), 0.889, 0.867, and 0.872 ( Table 2 ), respectively. Thus, the combined model integrating ultrasomics signatures and clinical features had a better performance in differentiation between HCC and ICC than the other two models. The combined model had more stable performance and higher generalization ability (p < 0.05).