The present research met with the approval of our Hospital Review Board (approval number KY2019-217) and written informed assent was exempted for all participators. The participants with HCC were enrolled from imaging database management system between March 2017 and March 2021. Gadolinium-ethoxybenzyl-diethylenetriamine-pentaacetic acid (Gd-EOB-DTPA, Primovist, Germany)-enhanced MRI imaging was carried out for a total of 284 patients who were suspected of HCC and experienced subsequent hepatic excision. The subject entry criteria were as follows: 1) Patients had a definite diagnosis of single primary HCC by pathology; 2) spanning less than 14 days between the first MRI examination and surgery. The exclusion criteria for patients included: 1) Patients had previous treatment, including partial hepatectomy, transcatheter arterial chemoembolization (TACE), radiotherapy, chemotherapy or needle biopsy (n = 23); 2) patients with tumor emboli in the bile duct or vessel (n = 24), or having distant transfer (n = 12); 3) Clinical and pathological data were incomplete (n = 10); 4) with poor image quality (n = 13). Finally, 202 patients (169 men and 33 women) were eligible and randomly divided into two datasets at a proportion of 7:3. Details of baseline, involving age, gender, tumor diameter, etiology of liver disease, cirrhosis, ascites, aspartate aminotransferase (AST), alanine aminotransferase (ALT), total bilirubin (TBIL), neutrophil count, lymphocyte count, neutrophil-to-lymphocyte ratio (NLR), albumin, serum α-fetoprotein (AFP), differentiation degree, were derived from medical records. Supplementary Method 1 presents the exhaustive MRI procedure. The flowchart for patient selection is demonstrated in Figure 1 .

All patients underwent subsequent radical resection of HCC. The detailed immunohistochemical protocol is explained in Supplementary Method 2 . Two observers without knowledge of the clinicopathologic parameters scored the slides independently using a previously validated scoring system. Any disagreements were arbitrated by a third rater. VEGF scoring was calculated using the immune reactive score (IRS) as depicted before (21–23): IRS = SI (staining intensity) × PP (percentage of positive cells). SI was graded as 0(negative), 1 (weak), 2 (moderate) and 3 (strong). PP was categorized as 0 (negative), 1 (less than 5% positive cells), 2 (6-25% positive cells), 3 (26-50% positive cells), 4 (51-75% positive cells), and 5 (more than 75% positive cells). The IRS was evaluated using ten visual fields from different regions of each tumor, and the average IRS was determined. The IRS classified VEGF status into four categories: grade 0 (IRS 0), grade 1 (IRS 1-4), grade 2 (IRS 5-9), and grade 3 (IRS 10-15). VEGF (grade 0 and grade 1) was defined as VEGF-negative expression, while VEGF (grade 2 and grade 3) was marked as VEGF-positive expression (24).

Two radiologists who were both lack of the histopathologic information analyzed the MR images on a Picture Archiving and Communication System (PACS) work-station. They negotiated every divergence to make a consensus. Supplemental Table 1 demonstrated the agreements of inter-observer about morphologic MRI features. The following qualitative morphologic variables were evaluated (20): (a) The tumor margin was classified as either smooth or irregular. It was defined irregular tumor margin if the ill-defined interface with normal tissue on HBP images was present; (b) arterial rim enhancement, described as a ring-enhanced with central low signal areas on arterial phase; (c) radiological capsule presence or absence.

All images of the HCC tumors were uploaded on the online Deepwise Research Platform (https://research.deepwise.com). Radiologist A, who was ignorant of the pathological conditions and had five years of expertise interpreting HCC imaging, delineated the HCC lesions. Radiologist B validated the region of interests (ROIs). ROIs were individually delineated on pictures from the arterial phase (AP), portal venous phase (PVP), balanced phase (BP), delayed phase (DP), and hepatobiliary phase (HBP). Twenty-two HCCs were arbitrarily chosen and re-segmented one week later by the two radiologists, in order to calculate the intra-class and inter-class correlation coefficient (ICCs).

We used logistic regression (univariate and multivariate) to evaluate the connection between basic information, morphologic features and VEGF expression in HCC of the training dataset. The stopping rule was a backward stepwise selection based on Akaike’s information criterion and the likelihood ratio test. In the multivariate analysis, the potential risk factors were selected as those with a p value < 0.05, and the clinical model derived from these variables was established for the training dataset.

The radiomics analysis technique involved ROI (HCCs) segmentation, feature extraction, feature selection, model building, and testing. The Python package Pyradiomics (Version 3.0.1, https://github.com/Radiomics/pyradiomics) was used for feature extraction (25). The following four steps presented the technique of feature extraction and selection in detail.

For additional radiomics signature development, logistic regression was used for radiomics characteristics determined by the F-test. For the optimal parameter selection in the training dataset, we utilized ten-fold cross-validation. The radiomics signature was created via the parameters that produced the best average area under the curve (AUC). The phases with AUCs larger than 0.700 were judged noteworthy and hence were chosen for the development of fusion radiomics signatures.

A likelihood ratio backward stepwise multivariable logistic regression model that combined clinical factors with radiomics signature was constructed to investigate whether these two variables enhanced prediction of VEGF status. Clinical variables with a significant correlation with VEGF, as well as radiomics signatures, were added into the multivariable combined model. The discriminatory ability of classifier model was analyzed by the AUCs, and the Delong test was used to compare AUCs. The Hosmer-Lemeshow test was used to determine the calibration performance between the combined model’s projected probability and the actual VEGF status.

Furthermore, we established a nomogram relying on a combined model that can simply calculate the likelihood of VEGF status. The clinical utility of the nomogram was appraised by means of decision curve analysis (28). By measuring the disparity between the true- and false-positive rate at different thresholds, net benefits were calculated (29). Figure 2 describes a diagram of a radiomics investigation.

We used Student’s t-test or Mann-Whitney U test for continuous variables and Chi-Square test or Fisher’s exact test for qualitative variables in training and test datasets. The statistical analyses were carried out using IBM SPSS Statistics 22.0 and MedCalc 20 software (MedCalc Software, Belgium). Statistical significance was defined as a two-tailed P-value < 0.05.

In Table 1 , demographics, laboratory results, and morphological characteristics of both the training and tested groups are summarized and compared. 202 individuals (169 male and 33 female) with solitary HCC were split into a training group (n = 142, 76 for VEGF positive, 66 for VEGF negative; 115 male and 27 female) and a testing group (n = 60, 32 for VEGF positive, 28 for VEGF negative; 54 male and 6 female).

Across the full cohort, there were 108 (53.5%) patients who were VEGF-positive and 94 (46.5%) patients who were VEGF-negative. VEGF status was not significantly different between the training and test datasets (p = 0.980). In the training cohort, univariate analysis revealed that serum ɑ-fetoprotein (AFP), irregular margin, and neutrophil-to-lymphocyte ratio (NLR) were significantly associated with VEGF (P < 0.05) ( Table 1 ). However, age, gender, tumor diameter, etiology of liver disease, cirrhosis, ascites, ALT, AST, TBIL, neutrophil count, lymphocyte count, albumin, differentiation degree, arterial rim enhancement, tumor capsule, and enhancement pattern were not statistically significant (P > 0.05). A backward stepwise multivariable logistic regression was used to examine every notable variable. Eventually, the meaningful variables for clinical model establishing were serum AFP level (odds ratio [OR] = 0.260; 95% confidence interval [CI]:0.107-0.633, P = 0.003), irregular tumor margin (p = 0.004, OR = 3.004, 95% CI: 1.434 - 6.295), and NLR (p = 0.027, OR = 0.622, 95% CI: 0.409 - 0.948) ( Table 2 ). In the training cohort, the clinical model’s AUC was 0.709 (95% CI: 0.624 - 0.794) while it was 0.725 (95% CI: 0.593 - 0.858) in the test cohort.

The features generated from arterial phase, PVP, balanced phase, delayed phase, and HBP were decreased to 1349, 1252, 1850, 1816, and 1202, respectively, after feature selection through the ICC test. Those robust features were applied for further analysis. Table 3 displays the forecast results for each of the five phases. Finally, the test set’s median 10-fold cross-validation AUC for PVP, HBP, AP, BP, and DP was 0.759, 0.731, 0.640, 0.592, and 0.692, respectively. The radiomics signatures of PVP and HBP (> 0.700) were captured to investigate further. Figure S1 presents the stability of radiomics characteristics in two phases (PVP and HBP) images. After F-test modeling for the prediction of VEGF-positive, 8 robust features in the PVP and 5 robust features in the HBP were regarded as the optimal feature set. The selected features were then used to create a logistic regression classifier model. Tables S2 , S3 provide more detailed information on the features. In the training cohort, the radiomics signature in the PVP obtained an AUC of 0.809 (95% CI: 0.740 - 0.878) while in the test dataset, it had an AUC of 0.759 (95% CI: 0.630 - 0.888). Likewise, in the training and test groups, the radiomics signature of HBP obtained an AUC of 0.792 (95% CI: 0.716 - 0.868) and 0.731 (95% CI: 0.597 - 0.865), respectively.

Furthermore, we propose a multiphase fusion model that combines the radiomics signatures (PVP and HBP) with a multivariable logistic regression model to study the maximum potential utilization of both radiomics signatures in diverse MRI phases. The fusion radiomics signature performed well in training and test groups, with AUCs of 0.892 (95% CI: 0.839 - 0.945) and 0.800 (95% CI: 0.682 - 0.918), respectively. When compared to the clinical model, the difference was significant in the training (p = 0.0003) cohort; however, the p value in the test dataset was 0.3759. We also noticed substantial variations in fusion signatures between VEGF (+) and VEGF (-) patients, and the p values of both (training and test) groups were less than 0.001.

The combined model performed best for VEGF prediction, combining all potential parameters including the fusion radiomics signature and major clinical characteristics (irregular tumor margin and serum AFP levels). Compared with the clinical model, the combined model represents superior performance in training [0.936(0.898 - 0.974) vs. 0.709(0.624-0.794); p < 0.0001] and test [0.836(0.728-0.944) vs. 0.725(0.593-0.858); p = 0.093] cohorts, albeit the difference was not statistically meaningful (p = 0.093) in test datasets. Furthermore, the accuracies of the combined model for forecasting the VEGF were 86.6% (training group) and 71.7% (test group) ( Table 3 ).

Comparisons among AUCs between every two model were showed in Table S4 . The stratified analysis also revealed strong results in both age and gender categories ( Table S5 ). Figure 3 displays the three ROC curves for different models. Figure 4 depicts the created nomogram. The calibration curves are displayed in Figure 5 , and there is high conformity in the two groups (training: P = 0.820; test: P = 0.508). Figure 6 depicts the findings of the decision curve analysis (DCA). The DCA for combined model indicated that when the probability of VEGF-positive was between 0.060-0.976 and 0.250-0.823 in the training and test datasets, respectively, using the nomogram to evaluate pathological VEGF-positive had added benefit compared with treating all patients as positive or negative. The coefficients and relative weights of 13 selected radiomics features of logistic regression in combined model are displayed on Table S6 . A correlation heat map of the selected 13 radiomics features in the combined model is depicted on Figure S2 .